# Quark production and thermalization of the quark-gluon plasma

Sergio Barrera Cabodevila<EMAIL_ADDRESS>Instituto Galego
de Física de Altas Enerxías IGFAE, Universidade de Santiago de Compostela,
E-15782 Galicia-Spain  Carlos A. Salgado<EMAIL_ADDRESS>Instituto
Galego de Física de Altas Enerxías IGFAE, Universidade de Santiago de
Compostela, E-15782 Galicia-Spain  Axencia Galega de Innovación (GAIN), Xunta
de Galicia, Galicia-Spain  Bin Wu<EMAIL_ADDRESS>Instituto Galego de Física de
Altas Enerxías IGFAE, Universidade de Santiago de Compostela, E-15782 Galicia-
Spain

###### Abstract

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\leftrightarrow
2$ 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. Schlichting:2019abc ; Berges:2020fwq ; Gelis:2021zmx  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 Jalilian-Marian:1996mkd ; Kovchegov:1998bi ;
Mueller:1999fp .

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
Berges:2013eia ; Epelbaum:2013ekf . However, the classical field approximation
is known to break down before the system thermalizes as the classical theory
is, in essence, non-renormalizable Epelbaum:2014yja  and the results under
this approximation eventually become dependent of the ultra-violet (UV) cutoff
Berges:2013lsa ; Epelbaum:2014mfa . 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/\alpha_{s}$ with $\alpha_{s}$ the strong coupling Mueller:2002gd .
Although such a transition cannot be confirmed by detailed calculations at the
lowest order beyond classical diagrams Wu:2017rry ; Kovchegov:2017way ,
switching from classical field simulations to the Boltzmann equation has
become a common practice in the studies of pre-equilibrium dynamics
Schlichting:2019abc ; Berges:2020fwq ; Gelis:2021zmx .

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\leftrightarrow 2$ processes and the $1\leftrightarrow 2$ processes. Two
variations of the LO Boltzmann equation have been presented in the literature
Baier:2000sb ; Arnold:2002zm . In ref. Baier:2000sb , the diffusion
approximation Mueller:1999pi  has been utilized to the $2\leftrightarrow 2$
kernel. This variation is equivalent to that in ref. Arnold:2002zm , referred
to as Effective Kinetic Theory (EKT), under the leading logarithmic
approximation.

The $2\leftrightarrow 2$ collision kernel in diffusion approximation
Mueller:1999pi , 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 Blaizot:2013lga ;
Blaizot:2014jna . Including only the $2\leftrightarrow 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$ Blaizot:2014jna . 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\leftrightarrow 2$ processes
are taken into account Semikoz:1994zp ; Semikoz:1995rd ; Scardina:2014gxa ;
Xu:2014ega ; Epelbaum:2015vxa . When the number-changing processes are
considered, different theories could behave differently. In the spin-0 scalar
theory studied in ref. Lenkiewicz:2019glw , the BEC could form transiently
before it vanishes during the thermalization process. Unexpected from the
parametric estimates in ref. Blaizot:2011xf , 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\leftrightarrow 2$ processes are included Blaizot:2016iir ;
BarreraCabodevila:2022jhi . Moreover, more elaborated analyses of infrared
modes beyond the scope of the kinetic theory have showed no evidence for a BEC
either Kurkela:2012hp .

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 Baier:1996kr ; Zakharov:1996fv ;
Baier:1998kq . When the $1\to 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\to 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. Baier:2000sb  (and first confirmed
by detailed numerical simulations in ref. Kurkela:2015qoa ). In a spatially
homogeneous, non-expanding system, the bottom-up thermalization also emerges
as the last stage of thermalization in initially very dilute cases
Kurkela:2014tea ; BarreraCabodevila:2022jhi .

Allowing the production of quarks and antiquarks delays the thermalization
process Blaizot:2014jna ; Kurkela:2018oqw ; Kurkela:2018xxd ; Du:2020zqg ;
Du:2020dvp . 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\leftrightarrow 2$ processes are considered in the BEDA Blaizot:2014jna
. The thermalization process is still delayed when the $1\leftrightarrow 2$
processes are additionally taken into account as studied using EKT
Kurkela:2018xxd ; Du:2020dvp . Moreover, gluons are found to achieve kinetic
equilibrium among themselves before quarks and antiquarks approach thermal
equilibrium Kurkela:2018xxd . 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 Kurkela:2018oqw ;
Kurkela:2018xxd ; Du:2020zqg ; Du:2020dvp .

In this paper, we present, for the first time, a full set of BEDA with both
the $2\leftrightarrow 2$ and $1\leftrightarrow 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
top-down 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 and summarize our
main results.

In sec. II.1, we incorporate the $1\leftrightarrow 2$ kernel using the LPM
splitting rates Baier:1996kr ; Zakharov:1996fv ; Baier:1998kq ; Arnold:2008zu
into the BEDA while only the $2\leftrightarrow 2$ kernel for gluons and
$N_{f}$ flavors of massless quarks and antiquarks were included in ref.
Blaizot:2014jna . In this way, we assemble a full set of BEDA as an extension
of that for pure gluon systems in ref. Baier:2000sb . Its collision kernels
contain the following space-time dependent quantities, calculated from the
phase-space distribution functions of partons: the jet quenching parameter
$\hat{q}$, the effective temperature
$T_{*}\propto\hat{q}/(\alpha_{s}m_{D}^{2})$ with $m_{D}^{2}$ the screening
mass squared, and two more for each quark flavor $i$, denoted by
$\mathcal{I}^{i}_{c}$ and $\mathcal{\bar{I}}^{i}_{c}$ in eq. (II.1). In terms
of the latter two quantities, a definition of the effective net quark chemical
potential $\mu_{*}^{i}$ is given for flavor $i$ in eq. (11), which equals the
net quark chemical potential after thermal equilibration.

In sec. II.2, we streamline the BEDA to investigate spatially homogeneous
systems. 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\leftrightarrow 2$ kernel. This extends
previous discussions on the onset of a BEC in pure gluon systems
Blaizot:2011xf ; Kurkela:2012hp ; Blaizot:2013lga ; Blaizot:2014jna ;
Blaizot:2016iir ; BarreraCabodevila:2022jhi  to generic systems including
additional $N_{f}$ flavors of massless quarks and antiquarks. Similar to the
observations in refs. Blaizot:2016iir ; BarreraCabodevila:2022jhi , we find
that both the gluon and quark distributions in the limit $p\to 0$ are
identical to those in thermal equilibrium with the temperature and the net
quark chemical potential given by $T_{*}$ and $\mu_{*}^{i}$, respectively. As
a result, the gluon flux at $p=0$, which is proportional to
$(\lim\limits_{p\to 0}pf^{g}-T_{*})$ with $f^{g}$ the gluon distribution
Blaizot:2014jna , always vanishes and, hence, no gluons are deposited into a
BEC.

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. BarreraCabodevila:2022jhi . For this purpose,
the initial gluon distribution is chosen to be the same as those in
Blaizot:2013lga ; Blaizot:2014jna ; BarreraCabodevila:2022jhi  while the
distributions for $N_{f}=3$ flavors of quarks and antiquarks vanish initially,
as detailed in sec. II.3. In this case, the time evolution of the system is
characterized solely by $\hat{q}$ and $m_{D}^{2}$ or, equivalently, $T_{*}$
with $\mathcal{I}^{i}_{c}=\mathcal{\bar{I}}^{i}_{c}\propto m_{D}^{2}$ and
$\mu_{*}^{i}=0$. This subsection also covers the corresponding thermal
equilibrium states, and our definition of the equilibration time.

In sec. III.1, 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}\ll 1$. Exactly like the pure gluon case
Schlichting:2019abc ; BarreraCabodevila:2022jhi , the system thermalizes
through three stages for $N_{f}=3$. And we find that quark production does not
modify the parametric forms of $\hat{q}$ and $m_{D}^{2}$ during each of these
stages. In Stage 1, both $\hat{q}$ and $m_{D}^{2}$ 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_{*}\sim\,Q$ and $\mu_{*}^{i}=0$ up to $p\sim 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 $\hat{q}$ and
$m_{D}^{2}$, albeit increasing over time, are all parametrically smaller than
those from hard gluons. In Stage 2, $m_{D}^{2}$ receives dominant
contributions from soft gluons and keeps increasing over time while $\hat{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\sim
f_{0}^{1/4}Q$, soft quarks and antiquarks become more abundant than hard ones.
Later on, $T_{*}$ keeps cooling while the number density of quarks and
antiquarks $\sim(\hat{q}t)^{3/2}$ grows more rapidly than that of soft gluons
$\sim T_{*}\hat{q}t$. 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_{*}\sim f_{0}^{1/3}Q$. In Stage 3, both $\hat{q}$
and $m_{D}^{2}$ 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.
Baier:2000sb ; Kurkela:2014tea ; Kurkela:2015qoa ; Kurkela:2018oqw ;
Kurkela:2018xxd ; Du:2020zqg ; Du:2020dvp ; BarreraCabodevila:2022jhi .

In sec. IV.1, we carry out parametric estimates for initially over-populated
systems with $f_{0}\gg 1$. We find that the two stages preceding
thermalization in the pure gluon case Schlichting:2019abc ;
BarreraCabodevila:2022jhi  are not affected by the production of quarks and
antiquarks. And their contributions to $\hat{q}$ and $m_{D}^{2}$ 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}\ll 1$: $\hat{q}$ and $m_{D}^{2}$
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_{*}\sim 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 $\sim 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 Berges:2012ev ;
Schlichting:2012es ; Kurkela:2012hp  and kinetic theory Blaizot:2011xf ;
AbraaoYork:2014hbk ; Kurkela:2014tea ; Du:2020dvp ; BarreraCabodevila:2022jhi
. Full thermalization is established when all the quantities exit their
scaling regions and approach their corresponding thermal equilibrium values.

Complementary quantitative studies are provided in sec. III.2 for initially
under-populated systems with $f_{0}=10^{-2},10^{-4}$ and $10^{-6}$ and in sec.
IV.2 for initially over-populated systems with $f_{0}=1,10$ and $100$. One
universal feature for the cases with $f_{0}\geq 10^{-4}$ that could not
revealed by the above parametric estimates is as follows. Before achieving
full thermalization, gluons undergo top-down 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. These observations qualitatively agree with those using EKT in refs.
Kurkela:2018oqw ; Kurkela:2018xxd ; Du:2020zqg ; Du:2020dvp . 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$.

Besides, the $1\to 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 Barrera_Cabodevila_hBEDA_2023
.

##  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 Baier:2000sb ;
BarreraCabodevila:2022jhi .

###  II.1 The BEDA for gluons, quarks and antiquarks

The leading-order QCD Boltzmann equation includes not only the
$2\leftrightarrow 2$ processes but also the medium-induced $1\leftrightarrow
2$ splitting Baier:2000sb ; Arnold:2002zm

$\displaystyle\left(\partial_{t}+{\bm{v}}\cdot\nabla_{\bm{x}}\right)f^{a}=$
$\displaystyle C^{a}_{2\leftrightarrow 2}+C^{a}_{1\leftrightarrow 2},$ (1)

where $a=g,q^{i}$ and $\bar{q}^{i}$ respectively standing for gluons, quarks
and antiquarks of the $i^{\text{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

$\displaystyle s=s_{g}+s_{q}+s_{\bar{q}},$ (2a) $\displaystyle
n=2N_{c}\int\frac{d^{3}p}{(2\pi)^{3}}\bigg{[}2C_{F}f^{g}+\sum\limits_{i=1}^{N_{f}}(f^{q_{i}}+f^{\bar{q}_{i}})\bigg{]}\equiv
n_{g}+n_{q}+n_{\bar{q}},$ (2b)
$\displaystyle\epsilon=2N_{c}\int\frac{d^{3}p}{(2\pi)^{3}}p\left[2C_{F}f^{g}+\sum\limits_{i=1}^{N_{f}}(f^{q_{i}}+f^{\bar{q}_{i}})\right]\equiv\epsilon_{g}+\epsilon_{q}+\epsilon_{\bar{q}},$
(2c)

where the number of colors $N_{c}=3$, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$, and the
entropy densities of gluons, quarks and antiquarks are respectively defined as

$\displaystyle
s_{g}\equiv-4N_{c}C_{F}\int\frac{d^{3}p}{(2\pi)^{3}}\left[f^{g}\log
f^{g}-(1+f^{g})\log(1+f^{g})\right],$ (3) $\displaystyle
s_{q}\equiv-2N_{c}\sum\limits_{i=1}^{N_{f}}\int\frac{d^{3}p}{(2\pi)^{3}}\left[f^{q^{i}}\log
f^{q^{i}}+(1-f^{q^{i}})\log(1-f^{q^{i}})\right],\qquad
s_{\bar{q}}=s_{q}\bigg{|}_{f^{q^{i}}\to f^{\bar{q}^{i}}}\,.$ (4)

In the diffusion approximation the $2\leftrightarrow 2$ kernel takes the form
Mueller:1999pi ; Hong:2010at ; Blaizot:2013lga ; Blaizot:2014jna

$\displaystyle C^{a}_{2\leftrightarrow
2}=\frac{1}{4}\hat{q}_{a}(t)\nabla_{{\bm{p}}}\cdot\left[\nabla_{{\bm{p}}}f^{a}+\frac{{\bm{v}}}{T^{*}(t)}f^{a}(1+\epsilon_{a}f^{a})\right]+\mathcal{S}_{a},$
(5)

where $\epsilon_{a}=1$ for bosons and $\epsilon_{a}=-1$ for fermions. Here,
the jet quenching parameter Baier:1996sk  for species $a$ is defined as
$\hat{q}_{a}=C_{a}\hat{\bar{q}}$, and

$\displaystyle\hat{\bar{q}}\equiv
4\pi\alpha_{s}^{2}\mathcal{L}\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\bigg{\\{}2N_{c}f^{g}\left(1+f^{g}\right)+\sum\limits_{i=1}^{N_{f}}[f^{q^{i}}(1-f^{q^{i}})+f^{\bar{q}^{i}}(1-f^{\bar{q}^{i}})]\bigg{\\}},$
(6)

where $C_{a}=N_{c}$ and $C_{F}$ respectively for $g$ and for
$q^{i}/\bar{q}^{i}$, $\mathcal{L}\equiv\ln\frac{\langle
p_{t}^{2}\rangle}{m_{D}^{2}}$, and the typical transverse momentum broadening
$\langle p_{t}^{2}\rangle=\frac{2}{3}\langle p^{2}\rangle$ for isotropic
systems. The effective temperature is given by Blaizot:2013lga

$\displaystyle
T_{*}(t)\equiv\frac{\hat{q}_{A}}{2\alpha_{s}N_{c}\mathcal{L}m_{D}^{2}},$ (7)

where the screening mass squared is defined as

$\displaystyle m_{D}^{2}=m_{D,g}^{2}+m_{D,q}^{2}+m_{D,\bar{q}}^{2}\equiv
4\pi\alpha_{s}\int\frac{d^{3}\bm{p}}{(2\pi)^{3}}\frac{1}{|\bm{p}|}\bigg{[}2N_{c}f^{g}+\sum\limits_{i=1}^{N_{f}}(f^{q^{i}}+f^{\bar{q}^{i}})\bigg{]}$
(8)

with $m_{D,a}^{2}$ standing for the contribution from parton $a$. Here,
$m_{D}^{2}$ is the same as the effective mass squared $m^{2}_{\text{eff,g}}$
in ref. Arnold:2002zm , which reduces to the thermal gluon mass squared or one
half of the Debye mass squared at $O(\alpha_{s})$ for a QGP in thermal
equilibrium Bellac:2011kqa .

By extending the source terms in ref. Blaizot:2014jna , one has, for each
quark flavor $i$,

$\displaystyle\mathcal{S}_{q^{i}}=\frac{2\pi\alpha_{s}^{2}C_{F}^{2}\mathcal{L}}{p}\bigg{[}\mathcal{I}^{i}_{c}f^{g}(1-f^{q^{i}})-\bar{\mathcal{I}}^{i}_{c}f^{q^{i}}(1+f^{g})\bigg{]},$
$\displaystyle\mathcal{S}_{\bar{q}^{i}}=\frac{2\pi\alpha_{s}^{2}C_{F}^{2}\mathcal{L}}{p}\bigg{[}\bar{\mathcal{I}}^{i}_{c}f^{g}(1-f^{\bar{q}^{i}})-\mathcal{I}^{i}_{c}f^{\bar{q}^{i}}(1+f^{g})\bigg{]},$
$\displaystyle\mathcal{S}_{g}=-\frac{1}{2C_{F}}\sum\limits_{i=1}^{N_{f}}({\mathcal{S}_{q^{i}}+\mathcal{S}_{\bar{q}^{i}}}),$
(9)

where one needs to define two more parameters for each quark flavor:

$\displaystyle\mathcal{I}^{i}_{c}=\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{p}[f^{g}+f^{q^{i}}+f^{g}(f^{q^{i}}-f^{\bar{q}^{i}})],$
$\displaystyle\bar{\mathcal{I}}^{i}_{c}=\int\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{p}[f^{g}+f^{\bar{q}^{i}}+f^{g}(f^{\bar{q}^{i}}-f^{q^{i}})].$
(10)

In terms of these two quantities, we define an effective net quark chemical
potential for the $i^{\text{th}}$ flavor

$\displaystyle{\mu_{*}^{i}}={T_{*}}\ln\frac{\mathcal{I}_{c}^{i}}{\mathcal{\bar{I}}^{i}_{c}},$
(11)

which reduces to the corresponding net quark potential in thermal equilibrium,
denoted by $\mu^{i}$, as a result of the conservation of the net quark number.
Besides, the $2\leftrightarrow 2$ kernel also conserves the total energy and
the total parton number.111 Without the $1\leftrightarrow 2$ processes, one
may need to take into account the number of the BEC for initially over-
populated systems Blaizot:2014jna  in order to preserve this conservation.

The $1\leftrightarrow 2$ processes are described by

$\displaystyle C^{a}_{1\leftrightarrow 2}=$
$\displaystyle\int_{0}^{1}dx\sum\limits_{b,c}\bigg{[}\frac{1}{x^{3}}\frac{\nu_{c}}{\nu_{a}}C^{c}_{ab}({\bm{p}}/{x};{\bm{p}},{\bm{p}}(1-x)/x)-\frac{1}{2}C^{a}_{bc}(\bm{p};x{\bm{p}},(1-x){\bm{p}})\bigg{]},$
(12)

where $\nu_{a}$ is the number of spin times color degrees of freedom for
parton $a$: $\nu_{a}=2N_{c}$ for $q$ or $\bar{q}$ and $\nu_{a}=2(N_{c}^{2}-1)$
for $g$, and

$\displaystyle C^{a}_{bc}(\bm{p};x{\bm{p}},(1-x){\bm{p}})\equiv\frac{dI_{a\to
bc}(p)}{dxdt}\mathcal{F}^{a}_{bc}(\bm{p};x{\bm{p}},(1-x){\bm{p}})$ (13)

with

$\displaystyle\mathcal{F}^{a}_{bc}(\bm{p};{\bm{l}},{\bm{k}})=f^{a}_{\bm{p}}(1+\epsilon_{b}f^{b}_{\bm{l}})(1+\epsilon_{b}f^{c}_{\bm{k}})-f^{b}_{\bm{l}}f^{c}_{\bm{k}}(1+\epsilon_{a}f^{a}_{\bm{p}}).$
(14)

Here, $x$ is the momentum fraction carried by particle $b$ for the process
$a\to\,bc$, ${\bm{l}}\equiv x{\bf p}$ and ${\bm{k}}\equiv(1-x){\bf p}$.
Following ref. Baier:2000sb , the splitting rates in the deep LPM regime
Baier:1998kq ; Arnold:2008zu  (see also Appendix A) are used in this paper:

$\displaystyle\frac{dI_{g\to
gg}(p)}{dxdt}=\frac{\alpha_{s}N_{c}}{\pi}\frac{(1-x+x^{2})^{5/2}}{(x-x^{2})^{3/2}}\sqrt{\frac{\hat{q}_{A}}{p}},$
$\displaystyle\frac{dI_{g\to
q\bar{q}}(p)}{dxdt}=\frac{\alpha_{s}}{4\pi}[x^{2}+(1-x)^{2}]\bigg{[}\frac{\frac{C_{F}}{C_{A}}-x(1-x)}{x(1-x)}\bigg{]}^{\frac{1}{2}}\sqrt{\frac{\hat{q}_{A}}{p}},$
$\displaystyle\frac{dI_{q\to gq}(p)}{dxdt}=\frac{dI_{\bar{q}\to
g\bar{q}}(p)}{dxdt}=\frac{\alpha_{s}C_{F}}{2\pi}\frac{1+(1-x)^{2}}{x}\bigg{[}\frac{1-x+\frac{C_{F}}{C_{A}}x^{2}}{x(1-x)}\bigg{]}^{\frac{1}{2}}\sqrt{\frac{\hat{q}_{A}}{p}},$
$\displaystyle\frac{dI_{q\to
qg}(p)}{dxdt}=\frac{dI_{\bar{q}\to\bar{q}g}(p)}{dxdt}=\frac{\alpha_{s}C_{F}}{2\pi}\frac{1+x^{2}}{1-x}\bigg{[}\frac{x+\frac{C_{F}}{C_{A}}(1-x)^{2}}{x(1-x)}\bigg{]}^{\frac{1}{2}}\sqrt{\frac{\hat{q}_{A}}{p}},$
(15)

where $q$ and $\bar{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

$\displaystyle C^{g}_{1\leftrightarrow 2}=$
$\displaystyle\int_{0}^{1}\frac{dx}{x^{3}}\bigg{[}C^{g}_{gg}+\frac{1}{2C_{F}}\sum\limits_{i=1}^{N_{f}}(C^{q^{i}}_{gq^{i}}+C^{\bar{q}^{i}}_{g\bar{q}^{i}})\bigg{]}({\bm{p}}/{x};{\bm{p}},{\bm{p}}(1-x)/x)$
$\displaystyle-\int_{0}^{1}dx\bigg{[}\frac{1}{2}C^{g}_{gg}+\sum\limits_{i=1}^{N_{f}}C^{g}_{q^{i}\bar{q}^{i}}\bigg{]}(\bm{p};x{\bm{p}},(1-x){\bm{p}}),$
(16) $\displaystyle C^{q^{i}}_{1\leftrightarrow 2}=$
$\displaystyle\int_{0}^{1}dx\bigg{[}\frac{1}{x^{3}}\bigg{(}C^{q^{i}}_{q^{i}g}+2C_{F}C^{g}_{q^{i}\bar{q}^{i}}\bigg{)}({\bm{p}}/{x};{\bm{p}},{\bm{p}}(1-x)/x)-C^{q^{i}}_{q^{i}g}(\bm{p};x{\bm{p}},(1-x){\bm{p}})\bigg{]},$
(17) $\displaystyle C^{\bar{q}^{i}}_{1\leftrightarrow 2}=$
$\displaystyle\int_{0}^{1}dx\bigg{[}\frac{1}{x^{3}}\bigg{(}C^{\bar{q}^{i}}_{\bar{q}^{i}g}+2C_{F}C^{g}_{\bar{q}^{i}q^{i}}\bigg{)}({\bm{p}}/{x};{\bm{p}},{\bm{p}}(1-x)/x)-C^{\bar{q}^{i}}_{\bar{q}^{i}g}(\bm{p};x{\bm{p}},(1-x){\bm{p}})\bigg{]}.$
(18)

It is straightforward to check that the $1\leftrightarrow 2$ kernel above
conserves the total energy and the net quark number for each flavor. As a
result, the BEDA with both the $2\leftrightarrow 2$ and $1\leftrightarrow 2$
kernels admits the following thermal fixed-point solutions:

$\displaystyle f^{g}_{eq}=\frac{1}{e^{\frac{p}{T}}-1},\qquad
f^{q^{i}}_{eq}=\frac{1}{e^{\frac{p-\mu^{i}}{T}}+1},\qquad
f^{\bar{q}^{i}}_{eq}=\frac{1}{e^{\frac{p+\mu^{i}}{T}}+1},$ (19)

where $T$ and $\mu^{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 $\mu^{i}$ can be obtained by the conservation of these two quantities.

###  II.2 The BEDA for spatially homogeneous systems

For spatially homogeneous systems, the BEDA reduces to

$\displaystyle\dot{f}^{a}=$
$\displaystyle\frac{1}{p^{2}}(p^{2}J^{a})^{\prime}+\mathcal{S}_{a}+C^{a}_{1\leftrightarrow
2},$ (20)

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

$\displaystyle\frac{1}{p^{2}}(p^{2}J^{a})^{\prime}\equiv\frac{1}{4}\frac{\hat{q}_{a}}{p^{2}}\bigg{[}p^{2}\bigg{(}{f^{a}}^{\prime}+\frac{1}{T_{*}}f^{a}(1+\epsilon_{a}f^{a})\bigg{)}\bigg{]}^{\prime}.$
(21)

As the $2\leftrightarrow 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

$\displaystyle\lim\limits_{p\to
0}p^{2}J^{a}=0=\lim\limits_{p\to\infty}p^{2}J^{a}.$ (22)

For gluons, the boundary condition $\lim\limits_{p\to 0}p^{2}J^{g}=0$ may not
be imposed consistently: If $f^{g}\propto 1/p$ at low $p$, it requires
$\lim\limits_{p\to 0}pf^{g}\to T_{*}$ Blaizot:2014jna . 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 Blaizot:2016iir ; BarreraCabodevila:2022jhi , 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\leftrightarrow 2$ kernel, evaluated with the splitting rates at leading
order in $x$ and $\mathcal{F}$ at leading order in $p$. To be more specific,
by assuming $p\ll p/x$ one has

$\displaystyle\mathcal{F}^{a}_{bc}\big{(}\frac{\bm{p}}{x};{\bm{p}},\frac{1-x}{x}{\bm{p}}\big{)}\approx
f^{a}_{\frac{\bm{p}}{x}}\big{(}1+\epsilon_{c}f^{c}_{\frac{\bm{p}}{x}}\big{)}+f^{b}_{\bm{p}}\big{(}\epsilon_{b}f^{a}_{\frac{\bm{p}}{x}}-f^{c}_{\frac{\bm{p}}{x}}\big{)}+pf^{b}_{\bm{p}}{f^{c}_{\frac{\bm{p}}{x}}}^{\prime}-\epsilon_{c}pf^{a}_{\frac{\bm{p}}{x}}{f^{c}_{\frac{\bm{p}}{x}}}^{\prime},$
(23)

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}\propto 1/p$. In this way, one has

$\displaystyle C^{g}_{1\leftrightarrow 2}$
$\displaystyle\approx\frac{\alpha_{s}}{2\pi}\sqrt{\frac{\hat{q}_{A}}{p}}\int\frac{dx}{x^{4}}\bigg{[}2N_{c}\mathcal{F}^{g}_{gg}+\sum\limits_{i=1}^{N_{f}}\big{(}\mathcal{F}^{q^{i}}_{gq^{i}}+\mathcal{F}^{\bar{q}^{i}}_{g\bar{q}^{i}}\big{)}\bigg{]}\bigg{(}\frac{\bm{p}}{x};{\bm{p}},\frac{1-x}{x}{\bm{p}}\bigg{)}$
$\displaystyle\approx\frac{\alpha_{s}N_{c}}{\pi}I_{a}\sqrt{\frac{\hat{q}_{A}}{p}}\frac{1}{p^{3}}\bigg{(}1-\frac{pf^{g}}{T_{*}}\bigg{)},$
$\displaystyle C^{q^{i}}_{1\leftrightarrow 2}$
$\displaystyle\approx\frac{\alpha_{s}C_{F}}{2\pi}\sqrt{\frac{\hat{q}_{F}}{p}}\int\frac{dx}{x^{3}}\bigg{[}\mathcal{F}^{q^{i}}_{q^{i}g}+\mathcal{F}^{g}_{q^{i}\bar{q}^{i}}\bigg{]}\bigg{(}\frac{\bm{p}}{x};{\bm{p}},\frac{1-x}{x}{\bm{p}}\bigg{)}$
$\displaystyle\approx{\alpha_{s}C_{F}}{\pi}(\mathcal{I}_{c}^{i}+\bar{\mathcal{I}}_{c}^{i})\sqrt{\frac{\hat{q}_{F}}{p}}\frac{1}{p^{2}}\left(\frac{1}{e^{-\frac{\mu^{i}_{*}}{T_{*}}}+1}-f^{q^{i}}\right),$
$\displaystyle C^{\bar{q}^{i}}_{1\leftrightarrow 2}$
$\displaystyle\approx{\alpha_{s}C_{F}}{\pi}(\mathcal{I}_{c}^{i}+\bar{\mathcal{I}}_{c}^{i})\sqrt{\frac{\hat{q}_{F}}{p}}\frac{1}{p^{2}}\left(\frac{1}{e^{\frac{\mu^{i}_{*}}{T_{*}}}+1}-f^{\bar{q}^{i}}\right)$
(24)

with $\mu_{*}^{i}$ defined in terms of $\mathcal{I}_{c}^{i}$ and
$\bar{\mathcal{I}}_{c}^{i}$ in eq. (11), and
$I_{a}\equiv\hat{q}_{A}/({\frac{4\alpha_{s}^{2}N_{c}^{2}}{\pi}\mathcal{L}})$.
Accordingly, in the limit $p\to 0$, the distribution functions take the form

$\displaystyle pf^{g}\to T_{*},\qquad
f^{q^{i}}\to\frac{1}{e^{-\frac{\mu^{i}_{*}}{T_{*}}}+1},\qquad
f^{\bar{q}^{i}}\to\frac{1}{e^{\frac{\mu^{i}_{*}}{T_{*}}}+1}.$ (25)

That is, they are the same as thermal distributions with the temperature and
the net quark chemical potential respectively given by $T_{*}$ and
$\mu_{*}^{i}$. As a result, one can consistently set $\lim\limits_{p\to
0}p^{2}J^{g}=0$, and there is no transient formation of a BEC according to
Blaizot:2014jna .

Neglecting the time dependence of the macroscopic quantities in eq. (II.2),
the BEDA admits the following approximate solutions

$\displaystyle pf^{g}$ $\displaystyle\approx
T_{*}\bigg{[}1-\bigg{(}1-\frac{p\left.f^{g}\right|_{t=0}}{T_{*}}\bigg{)}e^{-\left(\frac{p_{A*}}{p}\right)^{\frac{5}{2}}}\bigg{]},$
(26) $\displaystyle f^{q^{i}}$
$\displaystyle\approx\frac{1-e^{-\big{(}\frac{p_{F*}}{p}\big{)}^{\frac{5}{2}}}}{e^{-\frac{\mu^{i}_{*}}{T_{*}}}+1}+\left.f^{q^{i}}\right|_{t=0}e^{-\big{(}\frac{p_{F*}}{p}\big{)}^{\frac{5}{2}}},$
$\displaystyle f^{\bar{q}^{i}}$
$\displaystyle\approx\frac{1-e^{-\big{(}\frac{p_{F*}}{p}\big{)}^{\frac{5}{2}}}}{e^{\frac{\mu^{i}_{*}}{T_{*}}}+1}+\left.f^{\bar{q}^{i}}\right|_{t=0}e^{-\big{(}\frac{p_{F*}}{p}\big{)}^{\frac{5}{2}}}$

with

$\displaystyle p_{A*}\equiv(\hat{q}_{A}m_{D}^{4}t^{2}/4)^{\frac{1}{5}},\qquad
p_{F*}\equiv[\alpha_{s}C_{F}\pi(\mathcal{I}_{c}+\bar{\mathcal{I}}_{c})t]^{\frac{2}{5}}\hat{q}_{F}^{\frac{1}{5}}.$
(27)

###  II.3 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

$\displaystyle f=f^{g},\qquad F=f^{q},\qquad\bar{F}=f^{\bar{q}}$ (28)

with the quark flavor index neglected.

####  II.3.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. BarreraCabodevila:2022jhi , is affected by
allowing quark production, the initial gluon distribution is chosen to be the
same as that inspired by saturation physics in refs. Blaizot:2013lga ;
Blaizot:2014jna ; BarreraCabodevila:2022jhi :

$\displaystyle f(t=0)={f_{0}}\theta(Q-|\vec{p}|),\qquad
F(t=0)=F_{0}=0,\qquad\bar{F}(t=0)=\bar{F}_{0}=0.$ (29)

Correspondingly, the initial number density and energy density are
respectively given by

$\displaystyle
n(0)=\frac{N_{c}^{2}-1}{3\pi^{2}}{f_{0}}Q^{3},\quad\epsilon(0)=\frac{N_{c}^{2}-1}{4\pi^{2}}{f_{0}}Q^{4}.$
(30)

Since the BEDA is invariant under $F\leftrightarrow\bar{F}$, the system
preserves $F=\bar{F}$ in the subsequent evolution given $F=\bar{F}$ at the
initial time. Accordingly, one has $\mathcal{I}_{c}=\bar{\mathcal{I}}_{c}$ and
$\mu_{*}^{i}=0$. Furthermore, for $N_{f}=N_{c}=3$, one has

$\displaystyle
m_{D}^{2}=24\pi\alpha_{s}\mathcal{I}_{c}=24\pi\alpha_{s}\bar{\mathcal{I}}_{c}.$
(31)

That is, in this case the time evolution of the system is governed by the same
two quantities: the jet quenching parameter and $m_{D}^{2}$ as in pure gluon
systems.

For the initial distributions in eq. (29), the approximate solutions in eq.
(26) can be expressed as

$\displaystyle pf_{p}(t)$
$\displaystyle\approx\left\\{\begin{array}[]{ll}T_{*}&\text{for $p\lesssim
p_{A*}$}\\\
{pf_{0}\theta(Q-p)+[T_{*}-pf_{0}\theta(Q-p)]\left(\frac{p_{A*}}{p}\right)^{\frac{5}{2}}}&\text{for
$p\gtrsim p_{A*}$}\end{array}\right.,$ (34) $\displaystyle F_{p}(t)$
$\displaystyle=\bar{F}_{p}(t)=\frac{1}{2}[1-e^{-\left(\frac{p_{F*}}{p}\right)^{\frac{5}{2}}}]\approx\left\\{\begin{array}[]{ll}\frac{1}{2}&\text{for
$p\lesssim p_{F*}$}\\\
\frac{1}{2}\left(\frac{p_{F*}}{p}\right)^{\frac{5}{2}}&\text{for $p\gtrsim
p_{F*}$}\end{array}\right.$ (37)

with

$\displaystyle
p_{F*}\equiv(4\hat{q}_{A}m_{D}^{4}t^{2}/729)^{\frac{1}{5}}=({16}/{729})^{\frac{1}{5}}p_{A*}\approx
0.466p_{A*}.$ (38)

As $p_{F*}\sim p_{A*}$, we use

$\displaystyle p_{*}\equiv(\hat{q}_{A}m_{D}^{4}t^{2})^{\frac{1}{5}}\sim
p_{F*}\sim p_{A*}$ (39)

for the parametric estimates in the following sections,

####  II.3.2 The thermal equilibrium states and the equilibration time

As $\mu^{i}=0$, one can straightforwardly use energy conservation:

$\displaystyle\epsilon(0)=\epsilon_{eq}\equiv\epsilon_{eq,g}+\epsilon_{eq,q}+\epsilon_{eq,\bar{q}}=\bigg{[}2(N_{c}^{2}-1)+\frac{7}{4}N_{c}N_{f}+\frac{7}{4}N_{c}N_{f}\bigg{]}\frac{\pi^{2}}{30}T^{4}$
(40)

to obtain the thermal equilibrium temperature

$\displaystyle
T=\left[\frac{15(N_{c}^{2}-1)f_{0}}{4(N_{c}^{2}-1)+7N_{c}N_{f}}\right]^{\frac{1}{4}}\frac{Q}{\pi}\,=\left(\frac{60f_{0}}{16+\frac{21}{2}N_{f}}\right)^{\frac{1}{4}}\frac{Q}{\pi}\,.$
(41)

Then, by using the thermal distributions in eq. (19), one has

$\displaystyle
n_{eq}=n_{eq,g}+n_{eq,q}+n_{eq,\bar{q}}\equiv\big{[}2(N_{c}^{2}-1)+\frac{3}{2}N_{c}N_{f}+\frac{3}{2}N_{c}N_{f}\big{]}\frac{\zeta(3)}{\pi^{2}}T^{3}\,,$
$\displaystyle
s_{eq}=s_{eq,g}+s_{eq,q}+s_{eq,\bar{q}}\equiv\big{[}2(N_{c}^{2}-1)+\frac{7}{4}N_{c}N_{f}+\frac{7}{4}N_{c}N_{f}\big{]}\frac{2\pi^{2}}{45}T^{3},$
(42)

and

$\displaystyle\hat{q}_{A}=\frac{4\pi}{3}\alpha_{s}^{2}N_{c}\mathcal{L}\bigg{(}N_{c}+\frac{N_{f}}{2}\bigg{)}T^{3},\qquad
m_{D}^{2}=\frac{2\pi\alpha_{s}}{3}\bigg{(}N_{c}+\frac{N_{f}}{2}\bigg{)}T^{2}.$
(43)

And one can check that they satisfy
$\epsilon_{eq}=-P+Ts_{eq}+\sum\limits_{i=1}^{N_{f}}\mu^{i}(n_{eq,q^{i}}-n_{eq,\bar{q}^{i}})$
with $P=\epsilon_{eq}/3$. Note that for $N_{f}=3$ quarks and antiquarks
together account for over half of the amount of $\epsilon,n$ and $s$ while
they only contribute one third of the overall values of $\hat{q}_{A}$ and
$m_{D}^{2}$ 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 $\mathcal{M}$ as

$\big{|}1-\frac{\mathcal{M}(t_{eq})}{\mathcal{M}_{eq}}\big{|}=W\,,$ (44)

with $W$ a constant. Following Blaizot:2014jna , 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).

###  III.1 Parametric estimates for $f_{0}\ll 1$

In the pure gluon case, the system establishes thermal equilibrium through
three distinct stages for $f_{0}\ll 1$ Schlichting:2019abc ;
BarreraCabodevila:2022jhi . Below, we investigate whether and how the three
stages are modified when quark production is not ignored.

#### Stage 1. Overheating in the soft sector:
$0\ll\,Qt\ll\alpha_{s}^{-2}f_{0}$

During this stage, the properties of the system are dominated by hard gluons
with $p\sim Q$. And one has

$\displaystyle n_{g}\sim n_{h,g}\sim
f_{0}Q^{3},~{}~{}m_{D}^{2}\sim\alpha_{s}f_{0}Q^{2},~{}~{}\hat{q}\sim\alpha_{s}^{2}f_{0}Q^{3},~{}~{}T_{*}\sim
Q,$ (45)

where the number density of hard gluons is denoted by $n_{h,g}$ and
$\hat{q}_{A}$ is simply denoted by $\hat{q}$ here and below. Quarks,
antiquarks and soft gluons are predominantly generated by hard gluons through
the $1\leftrightarrow 2$ processes while, as justified below, the
$2\leftrightarrow 2$ processes play a parametrically negligible role in their
production.

Like the pure gluon case BarreraCabodevila:2022jhi , soft gluons are mostly
produced via the $g\leftrightarrow\;gg$ processes. From the $g\to gg$
splitting, the number density of gluons typically carrying momentum $p\ll Q$
can be estimated as

$\displaystyle
n_{s,g}(p)\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p}}n_{h,g}\sim\alpha_{s}^{2}f_{0}^{\frac{3}{2}}Q^{4}t\bigg{(}\frac{Q}{p}\bigg{)}^{\frac{1}{2}},$
(46)

and, equivalently,

$\displaystyle
f(p)\sim\frac{n_{s,g}(p)}{p^{3}}\sim\alpha_{s}^{2}f_{0}^{\frac{3}{2}}Qt\bigg{(}\frac{Q}{p}\bigg{)}^{\frac{7}{2}}.$
(47)

The above estimate for $f$ starts to break down when the reverse process:
$gg\to g$ becomes equally important. According to eq. (II.2), this occurs when
$pf(p)/T_{*}\sim 1$, corresponding to $p\sim
p_{*}=(\hat{q}_{A}m_{D}^{4}t^{2})^{\frac{1}{5}}\sim\alpha_{s}^{\frac{4}{5}}f_{0}^{\frac{3}{5}}(Qt)^{\frac{2}{5}}Q$.
At $p\lesssim p_{*}$, the $g\to\;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.
(34). Accordingly, the system witnesses a pronounced accumulation of soft
gluons typically carrying momentum $p_{s}\sim p_{*}$ with their number density
given by

$\displaystyle n_{s,g}\equiv
n_{s,g}(p_{s})\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p_{*}}}n_{h,g}\sim\alpha_{s}^{\frac{8}{5}}f_{0}^{\frac{6}{5}}(Qt)^{\frac{4}{5}}Q^{3}\sim
T_{*}p_{*}^{2}.$ (48)

The $2\leftrightarrow 2$ processes can be ignored in the above parametric
analysis: From multiple elastic scattering, the gluons typically pick up some
momentum broadening $\Delta
p\sim(\hat{q}t)^{\frac{1}{2}}\sim\alpha_{s}f_{0}^{\frac{1}{2}}(Qt)^{\frac{1}{2}}Q$,
which is parametrically negligible as $\Delta p\ll p_{*}$ at
$Qt\ll\alpha_{s}^{-2}f_{0}$.

Quarks and antiquarks can be produced via either the $g\to\;q/\bar{q}$
conversion in the $2\leftrightarrow 2$ scattering or the $g\to q\bar{q}$
process. According to the source terms in eq. (9), the $g\to\;q/\bar{q}$
conversion yields

$\displaystyle n_{q}^{g\to q}\sim
m_{D}^{4}t\sim\alpha_{s}^{2}f_{0}^{2}Q^{4}t.$ (49)

The quarks and antiquarks produced in this process typically carry hard
momentum $Q$ as hard gluons dominate $m_{D}^{2}$. From the $g\to\;q\bar{q}$
splitting, the number density of produced $q$ or $\bar{q}$ carrying a typical
momentum $p$ can be estimated as

$\displaystyle n^{g\to q\bar{q}}_{q}(p)\sim$ $\displaystyle
x\frac{d^{2}I_{g\to
q\bar{q}}(Q)}{dxdt}tn_{h,g}\sim\alpha_{s}\sqrt{x\frac{\hat{q}t^{2}}{Q}}n_{h,g}\sim\alpha_{s}^{2}f_{0}^{\frac{3}{2}}x^{\frac{1}{2}}Q^{4}t\qquad\text{with
$x\sim\frac{p}{Q}$},$ (50)

and, equivalently,

$\displaystyle F(p)\sim\frac{n^{g\to
q\bar{q}}_{q}(p)}{p^{3}}\sim\frac{1}{2}\left(\frac{p_{*}}{p}\right)^{\frac{5}{2}}.$
(51)

The above estimate using single branching starts to break down at $p\sim
p_{*}$ when $F(p_{*})\sim 1/(e^{-\frac{\mu^{i}_{*}}{T_{*}}}+1)=1/2$. For
$p\lesssim p_{*}$, the reverse process $q\bar{q}\to g$ becomes non-negligible
and its rivalry with the $g\to\;q\bar{q}$ process leads to thermalization of
this sector: $F\to 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\to q\bar{q}$ splitting are mostly hard with a number
density given by

$\displaystyle n_{q}^{g\to
q\bar{q}}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{Q}}n_{h,g}\ll
n_{s,g}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p_{*}}}n_{h,g}.$ (52)

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\sim p_{*}$ is the most
likely to be emitted via $g\to gg$. This difference is due to the soft
divergence of the latter process, which is absent in the former. As $f_{0}\ll
1$, the $g\to q\bar{q}$ splitting gives the dominant contribution, leading to
a linear growth of the quark/antiquark number over time:

$\displaystyle n_{q}\sim n_{q}^{g\to
q\bar{q}}\sim\alpha_{s}^{2}f_{0}^{\frac{3}{2}}Q^{4}t.$ (53)

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, $\hat{q}$ receives dominant contributions from hard gluons only.
For $m_{D}^{2}$, soft gluons and soft quarks with typical momentum $p_{s}\sim
p_{*}$ make the following contributions:

$\displaystyle
m_{Ds,g}^{2}\sim\alpha_{s}T_{*}p_{*}\sim\alpha_{s}^{\frac{9}{5}}f_{0}^{\frac{3}{5}}(Qt)^{\frac{2}{5}}Q^{2},\qquad
m_{Ds,q}^{2}\sim\alpha_{s}p_{*}^{2}\sim\alpha_{s}^{\frac{13}{5}}f_{0}^{\frac{6}{5}}(Qt)^{\frac{4}{5}}Q^{2}.$
(54)

They are both parametrically smaller than the contributions from hard gluons
with $m_{D}^{2}\sim\alpha_{s}f_{0}Q^{2}$ while the contribution from hard
quarks $\sim\alpha_{s}n_{q}/Q\sim\alpha_{s}^{3}f_{0}^{\frac{3}{2}}Q^{3}t$ is
parametrically even smaller.

#### Stage 2. Cooling and overcooling in the soft sector:
$\alpha_{s}^{-2}f_{0}\ll\,Qt\ll\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$

Figure 1: 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$.

During this stage, most of the partons in the system are still hard gluons,
leaving both $n_{g}\sim n_{h,g}\sim f_{0}Q^{3}$ and
$\hat{q}\sim\alpha_{s}^{2}n_{g}\sim\alpha_{s}^{2}f_{0}Q^{3}$ unmodified. On
the other hand, $m_{D}^{2}$ receives dominant contributions from soft gluons
whose momenta are pushed up to $p_{s}\sim(\hat{q}t)^{\frac{1}{2}}$ by multiple
elastic scattering. The number density of soft gluons can be estimated from
the $g\to gg$ splitting:

$\displaystyle
n_{s,g}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p_{s}}}n_{h,g}\sim
f_{0}^{\frac{1}{2}}(Qp_{s})^{\frac{3}{2}}\sim\alpha_{s}^{\frac{3}{2}}f_{0}^{\frac{5}{4}}(Qt)^{\frac{3}{4}}Q^{3}.$
(55)

And their contribution to $m_{D}^{2}$ takes the form

$\displaystyle
m_{D}^{2}\sim\alpha_{s}\frac{n_{s,g}}{p_{s}}\sim\hat{q}^{\frac{1}{2}}p_{s}^{\frac{1}{2}}\sim\alpha_{s}^{\frac{3}{2}}f_{0}^{\frac{3}{4}}(Qt)^{\frac{1}{4}}Q^{2},$
(56)

which is parametrically larger than that from hard gluons starting from
$Qt\sim\alpha_{s}^{-2}f_{0}$. Accordingly, one has

$\displaystyle
T_{*}\sim\hat{q}/(\alpha_{s}m_{D}^{2})\sim\alpha_{s}^{-\frac{1}{2}}f_{0}^{\frac{1}{4}}(Qt)^{-\frac{1}{4}}Q.$
(57)

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
BarreraCabodevila:2022jhi . In the range of $p$ dominated by the $g\to gg$
splitting, one has

$\displaystyle
f(p)\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p}}\frac{n_{h,g}}{p^{3}}\sim\frac{T_{*}}{p_{s}}\bigg{(}\frac{p_{s}}{p}\bigg{)}^{\frac{7}{2}}.$
(58)

It holds down to $p\sim p_{s}\sim(\hat{q}t)^{\frac{1}{2}}$, where
${p_{s}f(p_{s})}/{T_{*}}\sim 1$ and the reverse process $gg\to g$ becomes
equally important. Consequently, detailed balance is established in the soft
sector with $p\lesssim p_{s}$. At higher $p$, the estimate in eq. (58) breaks
down at $p\sim\tilde{p}$ with

$\displaystyle\tilde{p}=p_{s}\left(\frac{p_{s}}{T_{*}}f_{0}\right)^{-\frac{2}{7}}\sim\alpha_{s}^{\frac{4}{7}}f_{0}^{\frac{1}{7}}(Qt)^{\frac{2}{7}}Q,$
(59)

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\gtrsim
p\gtrsim\tilde{p}_{F}\equiv\alpha_{s}^{\frac{4}{3}}f_{0}^{\frac{1}{3}}(Qt)^{\frac{2}{3}}Q$,
quarks and antiquarks are mostly produced by hard gluons via the $g\to
q\bar{q}$ splitting with

$\displaystyle
F(p)\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{Q}}\bigg{(}\frac{p}{Q}\bigg{)}^{\frac{1}{2}}\frac{n_{h,g}}{p^{3}}\sim
f_{0}^{\frac{1}{2}}\bigg{(}\frac{Q}{p}\bigg{)}^{\frac{1}{2}}\bigg{(}\frac{p_{s}}{p}\bigg{)}^{2}\sim\alpha_{s}^{2}f_{0}^{\frac{3}{2}}\bigg{(}\frac{Q}{p}\bigg{)}^{\frac{5}{2}}Qt.$
(60)

For $p\lesssim\tilde{p}$, the intermediate portion of
$f\propto\;p^{-\frac{7}{2}}$ contributes to $F$ as well. However, its
contribution is parametrically negligible for $p\gtrsim\tilde{p}_{F}$. 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

$\displaystyle
F(p)\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p}}\frac{n_{g}(p)}{p^{3}}\sim\alpha_{s}^{2}\frac{\hat{q}t^{2}}{p}\frac{n_{h,g}}{p^{3}}\sim\frac{1}{2}\bigg{(}\frac{p_{s}}{p}\bigg{)}^{4}.$
(61)

By comparing this estimate with eq. (60), one can see that it becomes dominant
only for $p\lesssim\tilde{p}_{F}$ (note that one always has
$p_{s}\lesssim\tilde{p}_{F}\ll\tilde{p}$ throughout this stage). And the above
estimate based on single branching breaks down at $p\sim p_{s}$ as
$F(p_{s})\sim 1/2$. At $p\lesssim p_{s}$, the system maintains a balance
between the $g\to q\bar{q}$ process and its reverse process, corresponding to
$F\to 1/2$. Accordingly, unlike Stage 1, the system witnesses a pronounced
accumulation of soft quarks and antiquarks with $p\sim p_{s}$ besides hard
ones with $p\sim Q$:

$\displaystyle n_{s,q}^{1\leftrightarrow
2}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p_{s}}}n_{s,g}\sim p_{s}^{3},\qquad
n_{h,q}^{1\leftrightarrow
2}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{Q}}n_{h,g}\sim
f_{0}^{\frac{1}{2}}Qp_{s}^{2}.$ (62)

Moreover, starting at $Qt\sim\alpha_{s}^{-2}$ with $p_{s}\sim
f_{0}^{\frac{1}{2}}Q$, soft quarks and antiquarks with $p\sim p_{s}$ become
parametrically more abundant than hard ones.

Let us then look into quark production via the $2\leftrightarrow 2$ process.
From the conversion term in (9), one can obtain the quark distribution

$\displaystyle
F(p)\sim\alpha_{s}m_{D}^{2}t\frac{f_{p}}{p}\sim\alpha_{s}m_{D}^{2}t\frac{T_{*}}{p_{s}^{2}}\bigg{(}\frac{p_{s}}{p}\bigg{)}^{\frac{9}{2}}\sim\frac{1}{2}\bigg{(}\frac{p_{s}}{p}\bigg{)}^{\frac{9}{2}},$
(63)

which is negligible compared to the contribution from the $g\to q\bar{q}$
process for $p\gtrsim p_{s}$ in eq. (61). On the other hand, the quark number
produced via this process is given by

$\displaystyle n_{q}^{2\leftrightarrow 2}\sim m_{D}^{4}t\sim
p_{s}^{3}\sim\alpha_{s}^{3}f_{0}^{\frac{3}{2}}(Qt)^{\frac{3}{2}}Q^{3},$ (64)

which is parametrically the same as the contribution from the $g\to q\bar{q}$
process. Since the quark number density is parametrically smaller than
$n_{s,g}$ for $p_{s}\lesssim f_{0}^{\frac{1}{3}}Q$ at
$Qt\lesssim\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$, the quark and antiquark
production does not modify the parametric form of $\hat{q}$ or $m_{D}^{2}$.

In summary, even though soft gluons contribute dominantly to $m_{D}^{2}$ right
after $Qt\sim\alpha_{s}^{-2}f_{0}$, they do not contribute significantly to
quark production until $Qt\sim\alpha_{s}^{-2}$. Afterwards, the soft sector of
the system becomes overcooled with $T_{*}\lesssim T$, and soft quarks and
antiquarks are dominantly produced by soft gluons via both the $g\rightarrow
q\bar{q}$ and $gg\rightarrow q\bar{q}$ processes although there are less soft
gluons than hard ones in the system. As a result, one has

$\displaystyle
n_{q}\sim\left\\{\begin{array}[]{ll}\alpha_{s}^{2}f_{0}^{\frac{3}{2}}Q^{4}t&\text{at
$Qt\lesssim\alpha_{s}^{-2}$}\\\
\alpha_{s}^{3}f_{0}^{\frac{3}{2}}(Qt)^{\frac{3}{2}}Q^{3}&\text{at
$\alpha_{s}^{-2}\lesssim
Qt\lesssim\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$}\end{array}\right..$ (67)

#### Stage 3. Reheating and mini-jet quenching:
$\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}\ll\,Qt\ll\alpha_{s}^{-2}f_{0}^{-\frac{3}{8}}$

At $Qt\sim\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$, one has

$\displaystyle p_{s}\sim(\hat{q}t)^{\frac{1}{2}}\sim f_{0}^{\frac{1}{3}}Q\sim
T_{*},\qquad n_{q}\sim n_{s,g}\sim n_{h,g}\sim f_{0}Q^{3},\qquad
m_{D,q}^{2}\sim m_{D,g}^{2}\sim\alpha_{s}f_{0}^{\frac{2}{3}}Q^{2},$ (68)

which is consistent with a thermalized QGP of temperature $T_{*}\sim
f_{0}^{\frac{1}{3}}Q$. This is qualitatively different from the previous
stages as both $\hat{q}$ and $m_{D}^{2}$ are now determined by soft partons.
This can be further justified by the fact that the relaxation time
$t_{rel}\sim\alpha_{s}^{-2}T_{*}^{-1}$ 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_{\text{br}}\sim\alpha_{s}^{2}\hat{q}t^{2}$ Baier:2000sb ;
Baier:2001yt , which rapidly degrades into soft partons with $p_{s}\sim T_{*}$
Iancu:2015uja  via democratic branching Blaizot:2013hx  within a time of order
$t$. As the energy density of the QGP, sourced by hard partons, is given by
$T_{*}^{4}\sim\epsilon_{s}\sim p_{\text{br}}n_{h,g}$, one has

$\displaystyle T_{*}\sim\alpha_{s}^{4}f_{0}(Qt)^{2}Q,\qquad
m_{D}^{2}\sim\alpha_{s}T_{*}^{2}\sim\alpha_{s}^{9}f_{0}^{2}(Qt)^{4}Q^{2},\qquad\hat{q}\sim\alpha_{s}^{2}T_{*}^{3}\sim\alpha_{s}^{14}f_{0}^{3}(Qt)^{6}Q^{3}$
(69)

and

$\displaystyle n_{q}\sim n_{g}\sim\alpha_{s}^{12}f_{0}^{3}(Qt)^{6}Q^{3}.$ (70)

In the above parametric estimate, one can safely neglect hard quarks and
antiquarks as their number density
$n_{h,q}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{Q}}n_{h,g}\sim\alpha_{s}^{8}f_{0}^{\frac{5}{2}}(Qt)^{4}Q^{3}$
is parametrically smaller than $n_{h,g}\sim f_{0}Q^{3}$ for
$Qt\lesssim\alpha_{s}^{-2}f_{0}^{-\frac{3}{8}}$. And the soft sector remains a
thermalized QGP with a time-dependent temperature $T_{*}$ as the relaxation
time $t_{rel}\sim
1/(\alpha_{s}^{2}T_{*})\sim\alpha_{s}^{-6}f_{0}^{-1}(Qt)^{-2}Q^{-1}\lesssim t$
at $Qt\gtrsim\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$.

Eventually, the system thermalizes at
$Qt\sim\alpha_{s}^{-2}f_{0}^{-\frac{3}{8}}$ 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 Baier:2000sb  or without
Kurkela:2014tea ; BarreraCabodevila:2022jhi  longitudinal expansion. The only
distinction is that the soft thermal bath is now a QGP, composed of gluons as
well as quarks and antiquarks.

###  III.2 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. Following Blaizot:2014jna , $\mathcal{L}$ is kept
as a constant in our simulations. And we show the results for $\mathcal{L}=1$
and $\alpha_{s}=0.1$ below.

####  III.2.1 Initially dilute systems with $f_{0}=0.01$

Figure 2: The time evolution of different macroscopic quantities for
$f_{0}=0.01$. The left panel shows the results of $T_{*}$, $m_{D}^{2}$,
$\hat{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.

As a direct comparison with the pure gluon case studied in ref.
BarreraCabodevila:2022jhi , 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:
$(\alpha_{s}^{-2}f_{0},\alpha_{s}^{-2},\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}},\alpha_{s}^{-2}f_{0}^{-\frac{3}{8}})=(1,100,464,562)$.
Although there is no significant separation between the time scales after
$Qt\sim\alpha_{s}^{-2}$, the three stages: overheating, cooling/overcooling,
and reheating are still manifest in the time evolution of $T_{*}$ for
$N_{f}=0$ BarreraCabodevila:2022jhi  (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. 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 $\hat{q}_{A}$, $m_{D}^{2}$ and $T_{*}$ are all almost
indistinguishable from those in the pure gluon case before $Qt\sim 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\sim 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\sim 73$, similar to that during Stage
3 in the limit $f_{0}\ll 1$ BarreraCabodevila:2022jhi , while the late-time
evolution for $N_{f}=3$ becomes qualitatively different. At $Qt\sim 73$, one
has $T_{*}\approx 0.130Q$, which is still above its thermal equilibrium value
($T=0.107Q$). In the meanwhile, $\hat{q}$ reduces about 9% compared to that
for $N_{f}=0$ while $m_{D}^{2}$ is less sensitive to $N_{f}$, only dropping
about 3%. Such a trend continues later on, causing
$T_{*}\propto\hat{q}/(\alpha_{s}m_{D}^{2})$ 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)^{\frac{3}{2}}$, as shown in the right panel of fig. 2.

Figure 3: Contributions of gluons ($g$), quarks ($q$) and antiquarks
($\bar{q}$) to $\hat{q}$, $m_{D}^{2}$, $s$ and $n$ for $f_{0}=0.01$. Besides,
for $N_{f}=3$, the time derivatives of the number densities for $g$,
$q+\bar{q}$ and $g+q+\bar{q}$ via the $2\leftrightarrow 2$ processes and the
$1\leftrightarrow 2$ processes are presented in the lower center and lower
right panels, respectively. The lower right panel also shows
$\dot{n}_{g}+\dot{n}_{q}=\dot{n}_{g}+\dot{n}_{\bar{q}}$ via the
$1\leftrightarrow 2$ processes.

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\sim 73$,
the contributions from quarks and antiquarks to $\hat{q}$ (the upper left
panel) and $m_{D}^{2}$ (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
$\hat{q}$ and $m_{D}^{2}$ 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\leftrightarrow 2$ processes (the lower
right panel) and the $2\leftrightarrow 2$ processes (the lower center panel)
around this time.

Around $Qt\sim 87$, the gluon number density starts to decrease. Shortly
afterwards, the contributions from gluons to $\hat{q}$ and $m_{D}^{2}$ 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
$\hat{q}$ and $m_{D}^{2}$ develop a peak over time as gluons always give a
larger contribution than quarks 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. (II.3.2).

Figure 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.

The production of quarks and antiquarks is only partially responsible for the
rapid decrease in $\dot{n}_{g}$ via the $1\leftrightarrow 2$ processes shown
in the lower right panel of fig. 3. Another contributing factor to its
decrease is the increasing significance of the $gg\to g$ process, which
counteracts the $g\to gg$ process primarily responsible for the earlier
increase in $n_{g}$. Starting around $Qt=200$, one has
$\dot{n}_{g}^{1\leftrightarrow 2}+\dot{n}_{q}^{1\leftrightarrow
2}=\dot{n}_{g}^{1\leftrightarrow 2}+\dot{n}_{\bar{q}}^{1\leftrightarrow
2}\approx 0$. This implies that the excess of gluons is dominantly eliminated
by the $g\to q/\bar{q}$ conversion in the $2\leftrightarrow 2$ kernel and the
$g\to q\bar{q}$ splitting. As for the $g\to 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 have already reached their thermal
equilibrium values). This is consistent with the observation in the left panel
of fig. 4: the gluon distribution for $Qt\gtrsim 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\to 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. Kurkela:2018oqw .

Figure 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).

Let us now discuss thermalization more quantitatively in terms of the
equilibration time defined in eq. (44). For $N_{f}=0$, it takes
$t^{N_{f}=0}_{eq}\approx 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}\approx 424Q^{-1}$ for
$N_{f}=3$. However, the lower left panel of 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 $\epsilon_{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. (44). 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 Kurkela:2018oqw ;
Kurkela:2018xxd ; Du:2020zqg ; Du:2020dvp .

The lower right panel of fig. 5 focuses on the high-$p$ behavior of the gluon
distribution at $Qt\geq 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\leq 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 Maldacena:1997re ;
Witten:1998qj , where the UV modes equilibrate first Lin:2008rw ;
Balasubramanian:2010ce ; Balasubramanian:2011ur ; Wu:2012rib ; Wu:2013qi ;
Steineder:2013ana ; Keranen:2015mqc ; Attems:2018vuw ; Bernamonti:2018vmw .

####  III.2.2 Initially very dilute systems

The parametric estimates in sec. III.1 are valid only if $f_{0}\ll 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.

Figure 6: The time evolution of different macroscopic quantities for
$f_{0}=10^{-4}$. The left panel shows the results of $T_{*}$, $m_{D}^{2}$,
$\hat{q}$ and $n_{g}$ for $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.

Let us start with $f_{0}=10^{-4}$, which gives ($\alpha_{s}^{-2}f_{0}$,
$\alpha_{s}^{-2}$, $\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$,
$\alpha_{s}^{-2}f_{0}^{-\frac{3}{8}}$) = (0.01, 100, 2154, 3161). The detailed
time-evolution of $T_{*}$, $m_{D}^{2}$, $\hat{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 the soft
sector of the system experiences overheating and cooling. For $N_{f}=0$,
$T_{*}$ stops decreasing at $Qt\approx 196$, dropping about 16% below the
corresponding thermal equilibrium temperature. At this time, $\hat{q}$ and
$m_{D}^{2}$ for $N_{f}=3$ are respectively 4.5% and 1.4% below their
corresponding values for $N_{f}=0$. Up to this time, $m_{D}^{2}$ can be
approximated by the scaling law of $(Qt)^{1/4}$ for a wide range of $t$.
Meanwhile, following a linear growth, $n_{q}$ can be better approximated by
the scaling law of $(Qt)^{3/2}$ for a noticeable range of $t$.

Figure 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.

For $f_{0}=10^{-4}$, the late-time evolution process is still qualitatively
different for $N_{f}=0$ and $N_{f}=3$, as shown in the left panel of fig. 7.
After $Qt=196$, the pure gluon system continuously reheats up to the thermal
equilibrium temperature ($T=0.0443Q$), as expected in the parametric
estimates. For $N_{f}=3$, the ensuing thermalization pattern becomes more
complex. In this case, the first cooling phase stops at $Qt\approx 270$,
leaving $T_{*}$ 6.7% above the thermal equilibrium temperature ($T=0.0337Q$).
Then, following the trend in the pure gluon case, it gradually heats up to
$T_{*}=0.0372Q$ at $Qt=956$. Note that $n_{g}$ has already overshot its
thermal equilibrium value around this time, as shown by the ratio of
$n_{g}/n_{eq,g}$ in the left panel of fig. 7. Afterwards, the system undergoes
the second phase of cooling with $T_{*}$ decreasing toward its thermal
equilibrium value.

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$.

Figure 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).

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\to 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\leq 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.

Figure 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.

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
($\alpha_{s}^{-2}f_{0}$, $\alpha_{s}^{-2}$,
$\alpha_{s}^{-2}f_{0}^{-\frac{1}{3}}$, $\alpha_{s}^{-2}f_{0}^{-\frac{3}{8}}$)
= ($10^{-4}$, $10^{2}$, $10^{4}$, 17783). The late-time behavior of $T_{*}$ is
shown in the left panel of fig. 9.222We find that late-time results are robust
to values of the infrared cutoff $p_{min}$. However, one needs to choose
properly the value of $p_{min}$ in order to obtain correct early-time results,
especially for $m_{D}^{2}$.  It shows that the system overcools irrespective
of the values of $N_{f}$ with $T_{*}$ reaching its minimum at $Qt\approx 809$
for $N_{f}=0$ and at $Qt\approx 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 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^{N_{f}=0}_{eq}=1.19\times 10^{4}Q^{-1}$ 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\times 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.

###  IV.1 Parametric estimates for $f_{0}\gg 1$

In the limit $f_{0}\gg 1$, one initially has

$\displaystyle n_{g}\sim n_{h,g}\sim
f_{0}Q^{3},~{}~{}m_{D}^{2}\sim\alpha_{s}f_{0}Q^{2},~{}~{}\hat{q}\sim\alpha_{s}^{2}f_{0}^{2}Q^{3},~{}~{}T_{*}\sim
f_{0}Q.$ (71)

Like the pure gluon case Kurkela:2014tea ; Schlichting:2019abc ;
BarreraCabodevila:2022jhi , the system goes through the following two stages
towards achieving thermal equilibrium.

#### Stage 1. Soft gluon radiation and overheating:
$0\ll\,Qt\ll(\alpha_{s}f_{0})^{-2}$

During this stage, $\hat{q}$ and $m_{D}^{2}$ 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
$T_{*}\sim f_{0}Q\gg T\sim f_{0}^{\frac{1}{4}}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\to gg$
splitting, the number density of gluons typically carrying momentum $p$ is
given by

$\displaystyle n_{g}(p)\sim
x\frac{dI(Q)}{dxdt}{n_{h,g}}(1+f_{0})t\sim\alpha_{s}\sqrt{\frac{\hat{q}}{p}}f_{0}n_{h,g}t\sim\alpha_{s}^{-1}\sqrt{\frac{\hat{q}}{p}}\hat{q}t\qquad\text{with
$x\sim p/Q$},$ (72)

and, accordingly, the gluon distribution takes the following parametric form

$\displaystyle
f(p)\sim\frac{n(p)}{p^{3}}\sim\alpha_{s}^{-1}\sqrt{\frac{\hat{q}}{p^{7}}}\hat{q}t.$
(73)

The above estimate is valid for $p_{s}\lesssim
p\lesssim\tilde{p}\equiv(\hat{q}^{3}t^{2}/\alpha_{s}^{2}f_{0}^{2})^{1/7}\sim(\alpha_{s}^{2}f_{0}^{2}Qt)^{2/7}Q$
with the typical soft momentum $p_{s}$ given by the balancing condition:

$\displaystyle
1\sim\frac{p_{s}f(p_{s})}{T_{*}}\sim\sqrt{\frac{\hat{q}}{p_{s}^{5}}}m_{D}^{2}t\Leftrightarrow
p_{s}\sim
p_{*}=(\hat{q}m_{D}^{4}t^{2})^{\frac{1}{5}}\sim(\alpha_{s}^{2}\,f_{0}^{2}\,Qt)^{\frac{2}{5}}Q.$
(74)

For $p\lesssim p_{s}\sim 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 $\tilde{p}$ given above).
Accordingly, besides hard gluons, there is a pronounced accumulation of soft
gluons carrying a typical momentum of order $p_{s}$:

$\displaystyle
n_{s,g}\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{p_{s}}}f_{0}n_{h,g}\sim
T_{*}p_{s}^{2}\sim\alpha_{s}^{\frac{8}{5}}f_{0}^{\frac{8}{5}}(Qt)^{\frac{4}{5}}n_{h,g}.$
(75)

Evidently, it is parametrically smaller than the hard gluon number density at
$Qt\lesssim(\alpha_{s}f_{0})^{-2}$.

Quarks and antiquarks are produced much more slowly than soft gluons. Via the
$g\to q/\bar{q}$ conversion in the $2\leftrightarrow 2$ processes, one has
$n_{q}^{g\to q}$$\sim m_{D}^{4}t$ $\sim\alpha_{s}^{2}f_{0}^{2}Q^{4}t$. It is
parametrically the same as that produced via $g\to q\bar{q}$:

$\displaystyle n_{q}\sim n_{q}^{g\leftrightarrow
q\bar{q}}\sim\alpha_{s}\sqrt{\frac{\hat{q}}{Q}}n_{h,g}t\sim\alpha_{s}^{2}f_{0}^{2}Q^{4}t.$
(76)

That is, the $1\leftrightarrow 2$ processes are not more efficient than the
$2\leftrightarrow 2$ processes in producing quarks and antiquarks. On the
other hand, for $p\gtrsim p_{s}$, the quark distribution is mostly determined
by the $g\to q\bar{q}$ splitting:

$\displaystyle
F(p)\sim\alpha_{s}\sqrt{\frac{\hat{q}t^{2}}{Q}}\bigg{(}\frac{p}{Q}\bigg{)}^{\frac{1}{2}}\frac{n_{h,g}}{p^{3}}\sim\sqrt{\frac{\hat{q}t^{2}}{Q}}\bigg{(}\frac{p}{Q}\bigg{)}^{\frac{1}{2}}\frac{m_{D}^{2}Q}{p^{3}}\sim\frac{1}{2}\left(\frac{p_{s}}{p}\right)^{\frac{5}{2}}$
(77)

with $F(p_{s})\sim 1/2$. And at lower $p$, quarks and antiquarks fill a
thermal distribution with $F\to 1/2$.

During this stage, the contributions from soft gluons, quarks and antiquarks
to $m_{D}^{2}$ and $\hat{q}$ are all parametrically smaller than those from
hard gluons. At the end of this stage at $Qt\sim(\alpha_{s}f_{0})^{-2}$, one
has $p_{s}\sim\tilde{p}\sim Q$, $n_{s,g}\sim n_{h,g}\sim f_{0}Q^{3}\gg
n_{q}\sim Q^{3}$ and

$\displaystyle\hat{q}^{g}\sim\alpha_{s}^{2}f_{0}^{2}Q^{3}\gg\hat{q}^{q}\sim\alpha_{s}^{2}p_{s}^{3}\sim\alpha_{s}Q^{3},\qquad
m^{2}_{D,g}\sim\alpha_{s}f_{0}Q^{2}\gg
m^{2}_{D,q}\sim\alpha_{s}p_{s}^{2}\sim\alpha_{s}Q^{2}$ (78)

with $\hat{q}^{a}$ 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:
$(\alpha_{s}f_{0})^{-2}\ll\,Qt\ll\alpha_{s}^{-\frac{7}{4}}(\alpha_{s}f_{0})^{-\frac{1}{4}}$

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 BarreraCabodevila:2022jhi

$\displaystyle
p^{2}\sim\hat{q}t\sim\alpha_{s}^{2}n_{g}^{2}t/p^{3},\qquad\epsilon\sim\,f_{0}Q^{4}\sim\,n_{g}p$
(79)

yields

$\displaystyle
n_{g}\sim\frac{(\alpha_{s}f_{0})^{\frac{5}{7}}}{\alpha_{s}}\frac{Q^{3}}{(Qt)^{\frac{1}{7}}},\qquad
p\sim(\alpha_{s}f_{0})^{\frac{2}{7}}(Qt)^{\frac{1}{7}}Q,\qquad\,f\sim\frac{(\alpha_{s}f_{0})^{-\frac{1}{7}}}{\alpha_{s}}\frac{1}{(Qt)^{\frac{4}{7}}},$
$\displaystyle\hat{q}\sim(\alpha_{s}f_{0})^{\frac{4}{7}}\frac{Q^{3}}{(Qt)^{\frac{5}{7}}},\qquad
m_{D}^{2}\sim(\alpha_{s}f_{0})^{\frac{3}{7}}\frac{Q^{2}}{(Qt)^{\frac{2}{7}}},\qquad
T_{*}\sim\frac{(\alpha_{s}f_{0})^{\frac{1}{7}}}{\alpha_{s}}\frac{Q}{(Qt)^{\frac{3}{7}}}\sim
pf.$ (80)

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 Blaizot:2011xf ; Kurkela:2012hp ;
BarreraCabodevila:2022jhi .

Allowing the production of quarks and antiquarks does not modify the scaling
laws derived above. As in Stage 1, the $2\leftrightarrow 2$ and
$1\leftrightarrow 2$ processes are equally efficient to produce quarks and
antiquarks, yielding

$\displaystyle n_{q}\sim
m_{D}^{4}t\sim\alpha_{s}\sqrt{\frac{\hat{q}}{p}}n_{g}t\sim(\alpha_{s}f_{0})^{\frac{6}{7}}(Qt)^{\frac{3}{7}}Q^{3},$
(81)

and

$\displaystyle F\sim\frac{n_{q}}{p^{3}}\sim\frac{1}{2}.$ (82)

That is, parametrically, the majority of partons in the system consistently
fill thermal distributions up to
$p\sim\hat{q}t\sim(\alpha_{s}f_{0})^{\frac{2}{7}}(Qt)^{\frac{1}{7}}Q$.
Accordingly, one has

$\displaystyle\hat{q}^{q}\sim\alpha_{s}^{2}n_{q}\sim\alpha_{s}^{2}p^{3},\qquad
m_{D,q}^{2}\sim\alpha_{s}\frac{n_{q}}{p}\sim\alpha_{s}p^{2}.$ (83)

At $Qt\lesssim\alpha_{s}^{-\frac{7}{4}}(\alpha_{s}f_{0})^{-\frac{1}{4}}$, the
typical momentum of partons $p$ is always parametrically lower than $T_{*}$.
As a result, $\hat{q}$ and $m_{D}^{2}$ do not receive significant
contributions from quarks and antiquarks, and the system is hence dominated by
gluons.

At $Qt\sim\alpha_{s}^{-\frac{7}{4}}(\alpha_{s}f_{0})^{-\frac{1}{4}}$, all the
above quantities approach their final equilibrium values and full
thermalization is established around this time.

###  IV.2 Quantitative studies

In this subsection, we carry out simulations for initially over-populated
systems with $f_{0}>0.308$ for the initial distributions given in eq. (29). We
take $\alpha_{s}=0.1$ for $f_{0}=1$ and $10$, and $\alpha_{s}=0.01$ for
$f_{0}=100$ with $\mathcal{L}$ fixed to 1.

####  IV.2.1 Initially dense system with $f_{0}=1$

Figure 10: The time evolution of different macroscopic quantities for
$f_{0}=1$. The left panel shows the results of $T_{*}$, $m_{D}^{2}$, $\hat{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 behavior in the limit
$f_{0}\gg 1$.

Figure 11: The gluon (left) and quark (right) distributions at different times
for $f_{0}=1$. At low momentum, $pf$ is confirmed numerically to equal to
$T_{*}$ (for $Qt\gtrsim 10^{-2}$) and the quark distribution $F$ shows the
trend to approach 1/2. Both gluon and quark distributions fit very well the
thermal distributions (thermal) at the end of the evolution.

The time evolution of $\hat{q}$, $m_{D}^{2}$, $T_{*}$, $n_{g}$ and $n_{q}$ for
$f_{0}=1$ is shown in the two panels of Fig. 10. At early times, $\hat{q}$,
$m_{D}^{2}$, $T_{*}$ and $n_{g}$, as displayed in the left panel, are
insensitive to $N_{f}$. Neither of them experience significant changes before
$Qt\sim 1$. Besides, the right 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 $\hat{q}$, $m_{D}^{2}$, $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\to 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, $\hat{q}$,
$m_{D}^{2}$, $T_{*}$ and $n_{g}$ do not fit the scaling laws derived in the
limit $f_{0}\gg 1$. Similarly, the growth of $n_{q}$ is shown to slow down,
but it does not exhibit the scaling behavior of $(Qt)^{\frac{3}{7}}$.

Figure 12: The top-down thermalization 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_{*}$.

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 $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$.

####  IV.2.2 Initially very dense systems

Figure 13:  The time evolution of different macroscopic quantities for
$f_{0}=10$. The left panel shows the results of $T_{*}$, $m_{D}^{2}$,
$\hat{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 behavior in the
limit $f_{0}\gg 1$.

The time evolution of $\hat{q}$, $m_{D}^{2}$, $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 negligible role at
early times. In terms of $\hat{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\approx 1.4$, $\hat{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_{D}^{2}$ 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_{D}^{2}$ do not demonstrate convincingly
the scaling behaviors in the cooling stage for $f_{0}\gg 1$.

Figure 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_{*}$.

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}\geq 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 right panel of fig. 14.
For example, at $Qt=25$, the relative difference between $f$ and the
corresponding thermal distribution with temperature $T_{*}$ is smaller than 5%
for all the values of $p\lesssim 2.0T_{*}$. Eventually, gluons, quarks and
antiquarks all equilibrate around $t_{eq}=130Q^{-1}$, which is about 11.6
times that in the pure gluon case as determined by $T_{*}$, $n$ and $s$.

Figure 15: The time evolution of different macroscopic quantities for
$f_{0}=100$ and $\alpha_{s}=0.01$. The left panel shows the results of
$T_{*}$, $m_{D}^{2}$, $\hat{q}$ and $n_{g}$ for $N_{f}=0$ (dotted) and
$N_{f}=3$ (solid). The right panel shows $n_{q}$ over time. In both panels,
the expected scaling behaviors of these quantities (dashed) are also shown.

Let us now study the even denser system with $f_{0}=100$ and
$\alpha_{s}=0.01$, where $\hat{q}$ and $m_{D}^{2}$ are known to exhibit the
scaling behaviors in eq. (IV.1) for $N_{f}=0$ BarreraCabodevila:2022jhi . As
shown in the left panel of fig. 15, $\hat{q}$, $m_{D}^{2}$, $n_{g}$ and
$T_{*}$ for $N_{f}=3$ first experience the slowly varying, overheating stage.
Then, they approach the expected scaling behaviors. Before exiting the scaling
regions, they are all barely distinguishable from those in the pure gluon
case. Regarding $n_{q}$ as plotted in the right panel of Fig. 15, it initially
grows linear with time, and then approaches the anticipated scaling behavior
of $(Qt)^{\frac{3}{7}}$.

Figure 16: Self-similar scaling behaviors of the parton distributions for
$f_{0}=100$ and $\alpha_{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. 15 follow the scaling laws. The right
panels compare the distributions at $Qt=0.3$ and the others rescaled according
to eq. (85).

In our parametric estimates, it suffices to make the approximation: $pf\sim
T_{*}$ and $F\sim 1/2$ up to the typical momentum
$\bar{p}\equiv(\alpha_{s}f_{0})^{\frac{2}{7}}(Qt)^{\frac{1}{7}}Q$ in Stage 2.
Keeping only the range of $p\lesssim\bar{p}$ yields

$\displaystyle{(Qt)^{-\frac{3}{7}}}D_{s}(p/\bar{p})\equiv pf/Q\sim
T_{*}\sim\frac{(\alpha_{s}f_{0})^{\frac{1}{7}}}{\alpha_{s}}{(Qt)^{-\frac{3}{7}}},\qquad
F_{s}(p/\bar{p})\equiv F\sim 1/2.$ (84)

And solving exactly the BEDA only fills in quantitative details that are not
parametrically more important than the above parametric form. Consistent with
eq. (84), keeping only the dominant terms in the collision kernels under the
assumption $f\gg 1$ yields the following self-similar solutions

$\displaystyle
f={(Qt)^{-\frac{3}{7}}}\frac{Q}{p}D_{s}(p/\bar{p})\equiv{(Qt)^{-\frac{4}{7}}}f_{s}(\tilde{p}),\qquad
F=F_{s}(\tilde{p})$ (85)

with $\tilde{p}\equiv(Qt)^{-\frac{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. Blaizot:2011xf ; Kurkela:2012hp ; Schlichting:2019abc . To justify the
self-similar form of $F$, let us assume the ansatz
$F=(Qt)^{\alpha}F_{s}((Qt)^{\beta}p/Q)$ and take the $1\leftrightarrow 2$
kernel as a start. Keeping only the dominant term from the $g\to q\bar{q}$
splitting, one has $C^{q}_{1\leftrightarrow
2}=(Qt)^{-1}\tilde{C}^{q}(\tilde{p})$ with $\tilde{C}^{q}$ only dependent of
$\tilde{p}$. Then, matching it to the time derivative of $F$ yields $\alpha=0$
and $\beta=-1/7$. One can easily see that including the $2\leftrightarrow 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.

Figure 17: The top-down thermalization of gluons for $f_{0}=100$ and
$\alpha_{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_{*}$.

After exiting the scaling phase, the pure gluon system equilibrates around
$t^{N_{f}=0}_{eq}\approx 634Q^{-1}$ 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}\geq
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.

###### Acknowledgements.

We thank Jean-Paul Blaizot and Xiaojian Du for useful discussions. This work
is supported by European Research Council project ERC-2018-ADG-835105
YoctoLHC; by Spanish Research State Agency under project PID2020-119632GB-
I00; and by Xunta de Galicia (Centro singular de investigación de Galicia
accreditation 2019-2022), by European Union ERDF. B.W. acknowledges the
support of the Ramón y Cajal program with the Grant No. RYC2021-032271-I and
the support of Xunta de Galicia under the ED431F 2023/10 project. S.B.C.
acknowledges the support of the Axudas de apoio á etapa predoutoral program
(Ref. ED481A 2022/279).

##  Appendix A The $1\to 2$ splitting rates in the deep LPM regime

The splitting rates for $a\to bc$ in the deep LPM regime can be read out from
ref. Arnold:2008zu :

$\displaystyle\frac{dI_{a\to
bc}(p)}{dxdt}=\frac{(2\pi)^{3}}{\nu_{a}p}\gamma_{a\leftrightarrow
bc}(p;xp,(1-x)p)$ (86)

where

$\displaystyle\gamma_{g\leftrightarrow gg}(p;xp,(1-x)p)$
$\displaystyle=\frac{\sqrt{2}d_{A}C_{A}\alpha_{s}}{(2\pi)^{4}}\,m_{D}^{2}\,\hat{\mu}_{\perp}^{2}(1,x,1{-}x;{\rm
A},{\rm A},{\rm A})\,\frac{1+x^{4}+(1-x)^{4}}{x^{2}(1-x)^{2}}\,,$ (87a)
$\displaystyle\gamma_{q\leftrightarrow gq}(p;xp,(1-x)p)$
$\displaystyle=\frac{\sqrt{2}d_{F}C_{F}\alpha_{s}}{(2\pi)^{4}}\,m_{D}^{2}\,\hat{\mu}_{\perp}^{2}(1,x,1{-}x;{\rm
F},{\rm A},{\rm F})\,\frac{1+(1-x)^{2}}{x^{2}(1-x)}\,,$ (87b)
$\displaystyle\gamma_{g\leftrightarrow q\bar{q}}(p;xp,(1-x)p)$
$\displaystyle=\frac{\sqrt{2}d_{F}C_{F}\alpha_{s}}{(2\pi)^{4}}\,m_{D}^{2}\,\hat{\mu}_{\perp}^{2}(1,x,1{-}x;{\rm
A},{\rm F},{\rm F})\,\frac{x^{2}+(1-x)^{2}}{x(1-x)}$ (87c)

with $d_{F}=N_{c}$ and $d_{A}=N_{c}^{2}-1$, and

$\displaystyle\hat{\mu}_{\perp}^{2}(x_{1},x_{2},x_{3};s_{1},s_{2},s_{3})=$
$\displaystyle\frac{gT}{m_{D}}\left[\frac{1}{\pi}\,x_{1}x_{2}x_{3}\,\frac{p}{T}\mathcal{L}\right]^{1/2}\Biggl{[}\tfrac{1}{2}(C_{s_{2}}+C_{s_{3}}-C_{s_{1}})x_{1}^{2}$
$\displaystyle+\tfrac{1}{2}(C_{s_{3}}+C_{s_{1}}-C_{s_{2}})x_{2}^{2}+\tfrac{1}{2}(C_{s_{1}}+C_{s_{2}}-C_{s_{3}})x_{3}^{2}\Biggr{]}^{1/2}.$
(88)

Eliminating $m_{D}^{2}$ and $T$ in terms of
$\hat{q}=2\alpha_{s}\mathcal{L}m_{D}^{2}C_{A}T$ in the above equations yields
the splitting rates in eq. (15). And one can easily check that the rates for
$g\to gg$ and $q\to gq$ agree with those in ref. Baier:1998kq  in the deep LPM
regime. Note that these splitting rates are derived for a static medium. In a
generic medium, they remain applicable when the formation time of the
splitting is shorter than the typical time scale during which $\hat{q}$
undergoes substantial changes (see Iancu:2018trm  for a detailed discussion in
longitudinally boost-invariant systems).

##  Appendix B Numerical method

In the spatially homogeneous case, eq. (1) reduces to

$\displaystyle\partial_{t}f^{a}=$ $\displaystyle C^{a}_{2\leftrightarrow
2}+C^{a}_{1\leftrightarrow 2}.$ (89)

In order to solve it numerically, one can use the usual Euler integration
method. That is, given $f^{a}$ at $t$, the distribution at the next time step
$t+\Delta t$ is given by

$\displaystyle f^{a}(t+\Delta t)=f^{a}(t)+\Delta t\\{C^{a}_{2\leftrightarrow
2}[f(t)]+C^{a}_{1\leftrightarrow 2}[f(t)]\\}.$ (90)

And the discretization of the $2\leftrightarrow 2$ and $1\leftrightarrow 2$
kernels in our code is as follows.

###  B.1 Discretization of the $2\leftrightarrow 2$ kernel

As shown in eq. (5), in diffusion approximation the $2\leftrightarrow 2$
kernel can be expressed as that in the Fokker-Planck equation plus some source
term written in eq. (9). Here, one only needs to calculate derivatives and a
few integrals including $\hat{q}$, $m_{D}^{2}$, $\mathcal{I}_{c}$ and
$\bar{\mathcal{I}}_{c}$, which can be easily computed with, e.g., the
trapezoidal rule. In order to compute the derivatives, we define two different
momentum grids: $p_{s}$ and $p$ with the $i$-th point
$p_{i}=(p_{s,i}+p_{s,i-1})/2$ and $i=0,1,\cdots N$. Note $p_{s,-1}=0$ is not
include in the grid $p_{s}$. The distribution function $f$ is defined over the
grid $p$ while $J$, over the grid $p_{s}$. The value of the function
distribution on the grid $p_{s}$, denoted by $f_{s}$, is computed via a linear
interpolation of $f$ on $p$. For example, this definition allows one to
compute the first order derivative at the $i$-th point of the $p_{s}$ grid as

$\displaystyle\frac{\partial f}{\partial
p}\bigg{|}_{p_{s,i}}=\frac{f_{i+1}-f_{i}}{p_{i+1}-p_{i}}$ (91)

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
$\mathcal{S}$ on $p$, one has

$\displaystyle C^{a}_{2\leftrightarrow
2,i}=\frac{p_{s,i}^{2}J_{s,i}-p_{s,i-1}^{2}J_{s,i-1}}{\Delta
V_{i}}+\mathcal{S}_{i}$ (92)

with $\Delta V_{i}=\frac{1}{3}(p_{s,i}^{3}-p_{s,i-1}^{3})$.

One also needs to specify two boundary conditions:

$\displaystyle p_{s,-1}^{2}J_{s,-1}=0,\qquad p_{s,N}^{2}J_{s,N}=0$ (93)

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=\infty$.

###  B.2 Discretization of the $1\leftrightarrow 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

$\displaystyle C^{1\leftrightarrow
2}(p_{i})=\int_{0}^{1}dx\mathcal{K}(p_{i},x),$ (94)

where $\mathcal{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/\bar{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\leftrightarrow 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 than the CPU code for the same set
of initial conditions and parameters.

## References

  * (1) S. Schlichting and D. Teaney, _The First fm/c of Heavy-Ion Collisions_ , _Ann. Rev. Nucl. Part. Sci._ 69 (2019) 447 [1908.02113].
  * (2) J. Berges, M.P. Heller, A. Mazeliauskas and R. Venugopalan, _QCD thermalization: Ab initio approaches and interdisciplinary connections_ , _Rev. Mod. Phys._ 93 (2021) 035003 [2005.12299].
  * (3) F. Gelis, _Some Aspects of the Theory of Heavy Ion Collisions_ , _Rept. Prog. Phys._ 84 (2021) 056301 [2102.07604].
  * (4) J. Jalilian-Marian, A. Kovner, L.D. McLerran and H. Weigert, _The Intrinsic glue distribution at very small x_ , _Phys. Rev. D_ 55 (1997) 5414 [hep-ph/9606337].
  * (5) Y.V. Kovchegov and A.H. Mueller, _Gluon production in current nucleus and nucleon - nucleus collisions in a quasiclassical approximation_ , _Nucl. Phys. B_ 529 (1998) 451 [hep-ph/9802440].
  * (6) A.H. Mueller, _Toward equilibration in the early stages after a high-energy heavy ion collision_ , _Nucl. Phys. B_ 572 (2000) 227 [hep-ph/9906322].
  * (7) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Turbulent thermalization process in heavy-ion collisions at ultrarelativistic energies_ , _Phys. Rev. D_ 89 (2014) 074011 [1303.5650].
  * (8) T. Epelbaum and F. Gelis, _Pressure isotropization in high energy heavy ion collisions_ , _Phys. Rev. Lett._ 111 (2013) 232301 [1307.2214].
  * (9) T. Epelbaum, F. Gelis and B. Wu, _Nonrenormalizability of the classical statistical approximation_ , _Phys. Rev. D_ 90 (2014) 065029 [1402.0115].
  * (10) J. Berges, K. Boguslavski, S. Schlichting and R. Venugopalan, _Basin of attraction for turbulent thermalization and the range of validity of classical-statistical simulations_ , _JHEP_ 05 (2014) 054 [1312.5216].
  * (11) T. Epelbaum, F. Gelis, N. Tanji and B. Wu, _Properties of the Boltzmann equation in the classical approximation_ , _Phys. Rev. D_ 90 (2014) 125032 [1409.0701].
  * (12) A.H. Mueller and D.T. Son, _On the Equivalence between the Boltzmann equation and classical field theory at large occupation numbers_ , _Phys. Lett. B_ 582 (2004) 279 [hep-ph/0212198].
  * (13) B. Wu and Y.V. Kovchegov, _Time-dependent observables in heavy ion collisions. Part I. Setting up the formalism_ , _JHEP_ 03 (2018) 158 [1709.02866].
  * (14) Y.V. Kovchegov and B. Wu, _Time-dependent observables in heavy ion collisions. Part II. In search of pressure isotropization in the $\varphi^{4}$ theory_, _JHEP_ 03 (2018) 157 [1709.02868].
  * (15) R. Baier, A.H. Mueller, D. Schiff and D.T. Son, _’Bottom up’ thermalization in heavy ion collisions_ , _Phys. Lett. B_ 502 (2001) 51 [hep-ph/0009237].
  * (16) P.B. Arnold, G.D. Moore and L.G. Yaffe, _Effective kinetic theory for high temperature gauge theories_ , _JHEP_ 01 (2003) 030 [hep-ph/0209353].
  * (17) A.H. Mueller, _The Boltzmann equation for gluons at early times after a heavy ion collision_ , _Phys. Lett. B_ 475 (2000) 220 [hep-ph/9909388].
  * (18) J.-P. Blaizot, J. Liao and L. McLerran, _Gluon Transport Equation in the Small Angle Approximation and the Onset of Bose-Einstein Condensation_ , _Nucl. Phys. A_ 920 (2013) 58 [1305.2119].
  * (19) J.-P. Blaizot, B. Wu and L. Yan, _Quark production, Bose–Einstein condensates and thermalization of the quark–gluon plasma_ , _Nucl. Phys. A_ 930 (2014) 139 [1402.5049].
  * (20) D.V. Semikoz and I.I. Tkachev, _Kinetics of Bose condensation_ , _Phys. Rev. Lett._ 74 (1995) 3093 [hep-ph/9409202].
  * (21) D.V. Semikoz and I.I. Tkachev, _Condensation of bosons in kinetic regime_ , _Phys. Rev. D_ 55 (1997) 489 [hep-ph/9507306].
  * (22) F. Scardina, D. Perricone, S. Plumari, M. Ruggieri and V. Greco, _Relativistic Boltzmann transport approach with Bose-Einstein statistics and the onset of gluon condensation_ , _Phys. Rev. C_ 90 (2014) 054904 [1408.1313].
  * (23) Z. Xu, K. Zhou, P. Zhuang and C. Greiner, _Thermalization of gluons with Bose-Einstein condensation_ , _Phys. Rev. Lett._ 114 (2015) 182301 [1410.5616].
  * (24) T. Epelbaum, F. Gelis, S. Jeon, G. Moore and B. Wu, _Kinetic theory of a longitudinally expanding system of scalar particles_ , _JHEP_ 09 (2015) 117 [1506.05580].
  * (25) R. Lenkiewicz, A. Meistrenko, H. van Hees, K. Zhou, Z. Xu and C. Greiner, _Kinetic approach to a relativistic BEC with inelastic processes_ , _Phys. Rev. D_ 100 (2019) 091501 [1906.12111].
  * (26) J.-P. Blaizot, F. Gelis, J.-F. Liao, L. McLerran and R. Venugopalan, _Bose–Einstein Condensation and Thermalization of the Quark Gluon Plasma_ , _Nucl. Phys. A_ 873 (2012) 68 [1107.5296].
  * (27) J.-P. Blaizot, J. Liao and Y. Mehtar-Tani, _The thermalization of soft modes in non-expanding isotropic quark gluon plasmas_ , _Nucl. Phys. A_ 961 (2017) 37 [1609.02580].
  * (28) S. Barrera Cabodevila, C.A. Salgado and B. Wu, _Thermalization of gluons in spatially homogeneous systems_ , _Phys. Lett. B_ 834 (2022) 137491 [2206.12376].
  * (29) A. Kurkela and G.D. Moore, _UV Cascade in Classical Yang-Mills Theory_ , _Phys. Rev. D_ 86 (2012) 056008 [1207.1663].
  * (30) R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne and D. Schiff, _Radiative energy loss of high-energy quarks and gluons in a finite volume quark - gluon plasma_ , _Nucl. Phys. B_ 483 (1997) 291 [hep-ph/9607355].
  * (31) B.G. Zakharov, _Fully quantum treatment of the Landau-Pomeranchuk-Migdal effect in QED and QCD_ , _JETP Lett._ 63 (1996) 952 [hep-ph/9607440].
  * (32) R. Baier, Y.L. Dokshitzer, A.H. Mueller and D. Schiff, _Medium induced radiative energy loss: Equivalence between the BDMPS and Zakharov formalisms_ , _Nucl. Phys. B_ 531 (1998) 403 [hep-ph/9804212].
  * (33) A. Kurkela and Y. Zhu, _Isotropization and hydrodynamization in weakly coupled heavy-ion collisions_ , _Phys. Rev. Lett._ 115 (2015) 182301 [1506.06647].
  * (34) A. Kurkela and E. Lu, _Approach to Equilibrium in Weakly Coupled Non-Abelian Plasmas_ , _Phys. Rev. Lett._ 113 (2014) 182301 [1405.6318].
  * (35) A. Kurkela and A. Mazeliauskas, _Chemical equilibration in weakly coupled QCD_ , _Phys. Rev. D_ 99 (2019) 054018 [1811.03068].
  * (36) A. Kurkela and A. Mazeliauskas, _Chemical Equilibration in Hadronic Collisions_ , _Phys. Rev. Lett._ 122 (2019) 142301 [1811.03040].
  * (37) X. Du and S. Schlichting, _Equilibration of the Quark-Gluon Plasma at Finite Net-Baryon Density in QCD Kinetic Theory_ , _Phys. Rev. Lett._ 127 (2021) 122301 [2012.09068].
  * (38) X. Du and S. Schlichting, _Equilibration of weakly coupled QCD plasmas_ , _Phys. Rev. D_ 104 (2021) 054011 [2012.09079].
  * (39) P.B. Arnold and C. Dogan, _QCD Splitting/Joining Functions at Finite Temperature in the Deep LPM Regime_ , _Phys. Rev. D_ 78 (2008) 065008 [0804.3359].
  * (40) J. Berges, S. Schlichting and D. Sexty, _Over-populated gauge fields on the lattice_ , _Phys. Rev. D_ 86 (2012) 074006 [1203.4646].
  * (41) S. Schlichting, _Turbulent thermalization of weakly coupled non-abelian plasmas_ , _Phys. Rev. D_ 86 (2012) 065008 [1207.1450].
  * (42) M.C. Abraao York, A. Kurkela, E. Lu and G.D. Moore, _UV cascade in classical Yang-Mills theory via kinetic theory_ , _Phys. Rev. D_ 89 (2014) 074036 [1401.3751].
  * (43) S. Barrera Cabodevila, C.A. Salgado and B. Wu, “hBEDA.” https://github.com/cabodevila/hBEDA, Nov., 2023.
  * (44) J. Hong and D. Teaney, _Spectral densities for hot QCD plasmas in a leading log approximation_ , _Phys. Rev. C_ 82 (2010) 044908 [1003.0699].
  * (45) R. Baier, Y.L. Dokshitzer, A.H. Mueller, S. Peigne and D. Schiff, _Radiative energy loss and p(T) broadening of high-energy partons in nuclei_ , _Nucl. Phys. B_ 484 (1997) 265 [hep-ph/9608322].
  * (46) M.L. Bellac, _Thermal Field Theory_ , Cambridge Monographs on Mathematical Physics, Cambridge University Press (3, 2011), 10.1017/CBO9780511721700.
  * (47) R. Baier, Y.L. Dokshitzer, A.H. Mueller and D. Schiff, _Quenching of hadron spectra in media_ , _JHEP_ 09 (2001) 033 [hep-ph/0106347].
  * (48) E. Iancu and B. Wu, _Thermalization of mini-jets in a quark-gluon plasma_ , _JHEP_ 10 (2015) 155 [1506.07871].
  * (49) J.-P. Blaizot, E. Iancu and Y. Mehtar-Tani, _Medium-induced QCD cascade: democratic branching and wave turbulence_ , _Phys. Rev. Lett._ 111 (2013) 052001 [1301.6102].
  * (50) J.M. Maldacena, _The Large N limit of superconformal field theories and supergravity_ , _Adv. Theor. Math. Phys._ 2 (1998) 231 [hep-th/9711200].
  * (51) E. Witten, _Anti-de Sitter space and holography_ , _Adv. Theor. Math. Phys._ 2 (1998) 253 [hep-th/9802150].
  * (52) S. Lin and E. Shuryak, _Toward the AdS/CFT Gravity Dual for High Energy Collisions. 3. Gravitationally Collapsing Shell and Quasiequilibrium_ , _Phys. Rev. D_ 78 (2008) 125018 [0808.0910].
  * (53) V. Balasubramanian, A. Bernamonti, J. de Boer, N. Copland, B. Craps, E. Keski-Vakkuri et al., _Thermalization of Strongly Coupled Field Theories_ , _Phys. Rev. Lett._ 106 (2011) 191601 [1012.4753].
  * (54) V. Balasubramanian, A. Bernamonti, J. de Boer, N. Copland, B. Craps, E. Keski-Vakkuri et al., _Holographic Thermalization_ , _Phys. Rev. D_ 84 (2011) 026010 [1103.2683].
  * (55) B. Wu, _On holographic thermalization and gravitational collapse of massless scalar fields_ , _JHEP_ 10 (2012) 133 [1208.1393].
  * (56) B. Wu, _On holographic thermalization and gravitational collapse of tachyonic scalar fields_ , _JHEP_ 04 (2013) 044 [1301.3796].
  * (57) D. Steineder, S.A. Stricker and A. Vuorinen, _Probing the pattern of holographic thermalization with photons_ , _JHEP_ 07 (2013) 014 [1304.3404].
  * (58) V. Keranen and P. Kleinert, _Thermalization of Wightman functions in AdS/CFT and quasinormal modes_ , _Phys. Rev. D_ 94 (2016) 026010 [1511.08187].
  * (59) M. Attems, Y. Bea, J. Casalderrey-Solana, D. Mateos, D. Santos-Oliván, C.F. Sopuerta et al., _Paths to equilibrium in non-conformal collisions_ , _EPJ Web Conf._ 175 (2018) 07030.
  * (60) A. Bernamonti, F. Galli, R.C. Myers and J. Oppenheim, _Holographic second laws of black hole thermodynamics_ , _JHEP_ 07 (2018) 111 [1803.03633].
  * (61) E. Iancu, P. Taels and B. Wu, _Jet quenching parameter in an expanding QCD plasma_ , _Phys. Lett. B_ 786 (2018) 288 [1806.07177].