Quark production, Bose-Einstein condensates and thermalization of the quark-gluon plasma

TL;DR

This study models the thermalization process of quark-gluon plasma, revealing how initial gluon distributions influence Bose-Einstein condensate formation and showing quark production delays thermalization.

Contribution

It derives coupled transport equations for gluons and quarks, and explores the conditions leading to Bose-Einstein condensate formation during plasma thermalization.

Findings

Different thermalization patterns depending on initial gluon density.

Quark production extends the thermalization time significantly.

Bose-Einstein condensate formation depends on initial conditions.

Abstract

In this paper, we study the thermalization of gluons and N_f flavors of massless quarks and antiquarks in a spatially homogeneous system. First, two coupled transport equations for gluons and quarks (and antiquarks) are derived within the diffusion approximation of the Boltzmann equation, with only 2<-> 2 processes included in the collision term. Then, these transport equations are solved numerically in order to study the thermalization of the quark-gluon plasma. At initial time, we assume that no quarks or antiquarks are present and we choose the gluon distribution in the form f = f_0 theta (1-p/Q_s) with Q_s the saturation momentum and f_0 a constant. The subsequent evolution of systems may, or may not, lead to the formation of a (transient) Bose condensate, depending on the value of f_0. In fact, we observe, depending on the value of f_0, three different patterns: (a) thermalization without gluon Bose-Einstein condensate (BEC) for f_0 < f_{0t}, (b) thermalization with transient BEC for f_{0t} < f_0 < f_{0c}, and (c) thermalization with BEC for f_{0c} < f_0. The values of f_{0t} and f_{0c} depend on N_f. When f_0> 1 > f_{0c}, the onset of BEC occurs at a finite time t_c ~ 1/((alpha_s f_0)^2 Q_s). We also find that quark production slows down the thermalization process: the equilibration time for N_f = 3 is typically about 5 to 6 times longer than that for N_f = 0 at the same Q_s.

Tables

Table 1. Table 1: Squares of the 2 ↔ 2 ↔ 2 2 2\leftrightarrow 2 scattering amplitudes in QCD, with spins and colors of all four partons summed over. The dominant contributions of each process in diffusion approximation are given in the third column. The terms proportional to s 2 t 2 superscript 𝑠 2 superscript 𝑡 2 \frac{s^{2}}{t^{2}} or s 2 u 2 superscript 𝑠 2 superscript 𝑢 2 \frac{s^{2}}{u^{2}} contribute to the diffusion currents while the terms proportional to s t 𝑠 𝑡 \frac{s}{t} or s u 𝑠 𝑢 \frac{s}{u} only contribute to the source terms. Here, q 1 subscript 𝑞 1 q_{1} ( q ¯ 1 subscript ¯ 𝑞 1 \bar{q}_{1} ) and q 2 subscript 𝑞 2 q_{2} ( q ¯ 2 subscript ¯ 𝑞 2 \bar{q}_{2} ) represent quarks (antiquarks) of different flavors.

$ab\leftrightarrow cd$	${\left|{ℳ}_{cd}^{ab}\right|}^2/g^4$	In diffusion approximation
$q_1{q}_2\leftrightarrow q_1{q}_2$	$4N_c{C}_F\left(\frac{s^2+u^2}{t^2}\right)$	$8N_c{C}_F\frac{s^2}{t^2}$
$q_1{\overline{q}}_2\leftrightarrow q_1{\overline{q}}_2$
${\overline{q}}_1{q}_2\leftrightarrow {\overline{q}}_1{q}_2$
${\overline{q}}_1{\overline{q}}_2\leftrightarrow {\overline{q}}_1{\overline{q}}_2$
$q_1{q}_1\leftrightarrow q_1{q}_1$	$4N_c{C}_F\left(\frac{s^2+u^2}{t^2}+\frac{s^2+t^2}{u^2}\right)-8C_F\frac{s^2}{tu}$	$8N_c{C}_F\left(\frac{s^2}{t^2}+\frac{s^2}{u^2}\right)$
${\overline{q}}_1{\overline{q}}_1\leftrightarrow {\overline{q}}_1{\overline{q}}_1$
$q_1{\overline{q}}_1\leftrightarrow q_1{\overline{q}}_1$	$4N_c{C}_F\left(\frac{s^2+u^2}{t^2}+\frac{t^2+u^2}{s^2}\right)-8C_F\frac{u^2}{st}$	$8N_c{C}_F\frac{s^2}{t^2}$
$q_1{\overline{q}}_1\leftrightarrow q_2{\overline{q}}_2$	$4N_c{C}_F\left(\frac{t^2+u^2}{s^2}\right)$	0
$q_1{\overline{q}}_1\leftrightarrow gg$	$8N_c{C}_F^2\left(\frac{u}{t}+\frac{t}{u}\right)-8N_c^2{C}_F\frac{t^2+u^2}{s^2}$	$-8N_c{C}_F^2\left(\frac{s}{t}+\frac{s}{u}\right)$
$q_{1}g\leftrightarrow q_{1}g$	$-8N_c{C}_F^2\left(\frac{u}{s}+\frac{s}{u}\right)+8N_c^2{C}_F\frac{u^2+s^2}{t^2}$	$-8N_c{C}_F^2\frac{s}{u}+16N_c^2{C}_F\frac{s^2}{t^2}$
${\overline{q}}_{1}g\leftrightarrow {\overline{q}}_{1}g$
$gg\leftrightarrow gg$	$16N_c^2\left(N_c^2-1\right)\left(3-\frac{su}{t^2}-\frac{st}{u^2}-\frac{tu}{s^2}\right)$	$16N_c^2\left(N_c^2-1\right)\left(\frac{s^2}{t^2}+\frac{s^2}{u^2}\right)$

Equations

$$ϵ=3P=\left[16+\frac{21}{2}{N}_f\right]\frac{{\pi }^2}{30}{T}^4,$$

$$D_t{f}_{\mathbf{p}}^a\equiv \left(\frac{\partial }{\partial t}+\mathbf{v}\cdot {\nabla }_{\mathbf{x}}\right)f_{\mathbf{p}}^a=\mathcal{C}[f_{\mathbf{p}}^a],$$

$$\mathcal{C}[f_{\mathbf{p}}^a]=$$

$$\times \left[f_{\mathbf{k}}^c{f}_{{\mathbf{k}}^{\prime }}^d(1+{ϵ}_a{f}_{\mathbf{p}}^a)(1+{ϵ}_b{f}_{{\mathbf{p}}^{\prime }}^b)-f_{\mathbf{p}}^a{f}_{{\mathbf{p}}^{\prime }}^b(1+{ϵ}_c{f}_{\mathbf{k}}^c)(1+{ϵ}_d{f}_{{\mathbf{k}}^{\prime }}^d))\right],$$

$$D_{t}f=-{\nabla }_{\mathbf{p}}\cdot \mathcal{J},$$

$$\mathcal{L}\simeq {\int }_{q_{min}}^{q_{max}}\frac{dq}{q},$$

$$D_{t}f=$$

$$D_{t}F=$$

$${\mathcal{J}}_g$$

$${\mathcal{J}}_q$$

$${\mathcal{S}}_g=-\frac{N_f}{C_F}{\mathcal{S}}_q=\frac{4\pi {\alpha }_s^2{C}_F{N}_f\mathcal{L}{\mathcal{I}}_c}{p}\left[F(1+f)-f(1-F)\right],$$

$${\mathcal{I}}_a$$

$${\mathcal{I}}_b$$

$${\mathcal{I}}_c$$

$$|{\mathcal{M}}_{gq}^{gq}{|}^2/g^4=-8N_c{C}_F^2\left(\frac{u}{s}+\frac{s}{u}\right)+8N_c^2{C}_F\frac{u^2+s^2}{t^2},$$

$$\tau =\frac{2{\alpha }_s^2{N}_c\mathcal{L}}{\pi }t,$$

$$\dot{f}=-\frac{1}{p^2}{\left(p^2{J}_g\right)}^{\prime }+\frac{C_F{N}_f}{N_c}{S}_g=-\frac{1}{4\pi p^2}{\mathcal{F}}_g^{\prime }-\frac{C_F{N}_f}{N_c}{S}_q,$$

$$\dot{F}=-\frac{C_F}{N_c}\frac{1}{p^2}{\left(p^2{J}_q\right)}^{\prime }+\frac{C_F^2}{N_c}{S}_q=-\frac{C_F}{N_c}\frac{1}{4\pi p^2}{\mathcal{F}}_q^{\prime }+\frac{C_F^2}{N_c}{S}_q,$$

$$\frac{{\mathcal{F}}_g}{4\pi p^2}\equiv J_g\equiv -I_a{f}^{\prime }-I_{b}f(1+f),$$

$$\frac{{\mathcal{F}}_q}{4\pi p^2}\equiv J_q\equiv -I_a{F}^{\prime }-I_{b}F(1-F),$$

$$S_g=-S_q=\frac{I_c}{p}\left[F(1+f)-f(1-F)\right].$$

$$I_a=2{\pi }^2{\mathcal{I}}_a={\int }_0^{\infty }dp{p}^2\left[N_{c}f(1+f)+N_{f}F(1-F)\right],$$

$$I_b=2{\pi }^2{\mathcal{I}}_b=2{\int }_0^{\infty }dpp\left(N_{c}f+N_{f}F\right),I_c=2{\pi }^2{\mathcal{I}}_c={\int }_0^{\infty }dpp\left(f+F\right).$$

$$n=4N_c\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}\left(C_{F}f+N_{f}F\right)\equiv n_g+n_q,$$

$$ϵ=4N_c\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}p\left(C_{F}f+N_{f}F\right)\equiv {ϵ}_g+{ϵ}_q.$$

$$s_g\equiv -4N_c{C}_F\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}\left[f\log f-(1+f)\log (1+f)\right],$$

$$s_q\equiv -4N_c{N}_f\int \frac{d^3\mathbf{p}}{(2\pi {)}^3}\left[F\log F+(1-F)\log (1-F)\right],$$

$$s=s_g+s_q.$$

$$\dot{n}=$$

$$\dot{ϵ}=$$

$$\dot{s}=$$

$$s^{+}\equiv$$

$$f=\frac{c_{-1}}{p}-\frac{1}{2}+\cdots ,$$

$$F=\frac{1}{2}+\cdots ,$$

$${\mathcal{F}}_g{|}_{p=0}=4\pi c_{-1}(I_a-I_b{c}_{-1})\text{mbox}and{\mathcal{F}}_q{|}_{p=0}=0.$$

$$\dot{n}=$$

$$T^{\ast }\equiv \frac{I_a}{I_b},$$

Full text

I Introduction

Understanding how a dense system of gluons evolves into a thermalized quark-gluon plasma (QGP) is an important, and theoretically challenging, problem. After two colliding nuclei pass through each other in a relativistic heavy ion collision (HIC), a dense system of gluons is believed to be produced in a time scale of order $t\sim 1/{Q_{s}}$, with $Q_{s}$ the saturation momentum characterizing the initial nuclear wave functions Mueller:1999:initial . In this early stage, $f_{0}$, the occupation number of the produced gluons with $p\lesssim Q_{s}$, may be as large as $1/{\alpha_{s}}$, where $\alpha_{s}$ is the strong coupling constant. Under such conditions, it has been argued that a Bose-Einstein condensate (BEC) may develop during the approach to equilibrium, provided inelastic, number changing, processes do not play a too important role Blaizot:2011xf ; Blaizot:2013:BEC . The effect of such inelastic processes remains a somewhat controversial issue. Of course, number changing processes exclude the existence of a BEC in the equilibrium state. The real issue is therefore whether a transient BEC can emerge as the system evolves towards thermalization. Various arguments against this possibility are presented in Ref. Kurkela:2012hp , while the calculations in Ref. Huang:2013:2to3 suggest that inelastic processes could amplify the growth of soft gluon modes, thereby accelerating the formation of a BEC Huang:2013:2to3 . We shall not attempt to resolve this issue here, but consider rather the effect of another type of inelastic processes leading to the variation in the gluon number, namely processes that involve the creation of quark-antiquark pairs.

The partons that are produced in the early stage of HIC are mostly gluons: the number of quarks and antiquarks is initially negligible compared to the large number of gluons. However, in a thermalized quark-gluon plasma, the energy density is given by

[TABLE]

where we have assumed non-interacting quarks and gluons, $N_{f}$ flavors of massless quarks (and antiquarks), and $T$ is the temperature. At the energies of RHIC and LHC, one may take $N_{f}=3$. In this case quarks and antiquarks carry $66\%$ of the total energy density. Therefore, the study of quark production in a dense system of gluons is obviously of great importance to fully understand the thermalization of the quark-gluon plasma.

In this paper, we obtain two coupled kinetic equations for both gluons and quarks (and antiquarks), using the Boltzmann equation in the diffusion approximation PhysicalKinetics . The collision term contains all the $2\leftrightarrow 2$ scatterings between quarks and gluons, but only those $2\leftrightarrow 2$ scatterings, with the exclusion of, for instance, inelastic $2\leftrightarrow 3$ processes. We assume the dominance of small angle scatterings which justifies the diffusion approximation. The baryon number density is assumed to be zero. As a result quarks and antiquarks are described by the same transport equation, which is coupled to that for gluons. These transport equations are solved numerically to study the thermalization of the quark-gluon plasma.

The present study complements that carried out in Ref. Blaizot:2013:BEC where quark production was ignored. As in Blaizot:2013:BEC , the discussion relies on the Boltzmann equation in the small angle approximationMueller:1999:Boltzmann ; Venugopalan:2000:Thermalization ; Blaizot:2013:BEC , and both quarks and gluons are taken to be massless. As in Blaizot:2013:BEC , we restrict ourselves to the study of a spatially homogeneous non-expanding system. In contrast to Ref. Blaizot:2013:BEC , we are able to follow, albeit approximately, the evolution of the system across the onset of BEC all the way to thermalization. This is achieved by imposing a specific boundary condition on the solution of the coupled equations at zero momentum. It is shown in Ref. Blaizot:2013:BEC that the formation of BEC starts in an over-populated system at a finite time $t_{c}$ when the gluon distribution $f$ becomes singular at $p=0$. In this paper, we show that, for $t>t_{c}$, no solution of the transport equations exists if the total number of partons with $p>0$ is required to be conserved. However, we find solutions by properly imposing a boundary condition that corresponds to a non-vanishing gluon flux at $p=0$. Those solutions are used to describe the evolution of the system beyond $t_{c}$ all the way to thermal equilibrium, with the number density of condensed particles being deduced from the gluon flux at $p=0$. Note that the procedure just outlined represents presumably a crude approximation to the actual dynamics of particles in the presence of a condensate, but it has the virtue of allowing us to follow continuously the system all the way to its actual thermal equilibrium state.

Quark production decreases the total number of gluons in the system and could potentially hinder the formation of a BEC. However the $2\leftrightarrow 2$ processes included in the Boltzmann equation conserve the total number of partons. As a result of this conservation law, a chemical potential develops dynamically as the system evolves ***In fact, because the thermalization of quarks proceeds at a slower pace than that of the (soft) gluons, two different chemical potentials develop dynamically, one for the quarks and one for the gluons. These chemical potentials converge to a common value only close to thermalization. The equilibrium state is achieved for a negative value of this chemical potential, provided the initial number of gluons is not too large. We qualify this situation as under-population. If, on the contrary, the initial population of gluons is large enough, no equilibrium exists without a BEC: this is the situation of over-population, which was found to occur in the absence of quark production, and was thoroughly studied in Blaizot:2013:BEC . Thus the present study shows that quark production delays the onset of BEC but does not prevent the occurrence of the phenomenon. In fact, because the growth of the population of soft gluon modes is a fast phenomenon, and quark production is relatively slow, one even encounters situations where a transient BEC appears in the course of the evolution to equilibrium, before being eventually suppressed when quark production takes over and eliminates the excess gluons prior to thermalization.

The paper is organized as follows. The transport equations for quarks and gluons are derived in Sec. II. In Sec. III, the parameters that characterize the thermodynamic equilibrium are determined from the initial conditions, assuming that the total parton number is fixed. Our main results obtained by solving the transport equations for various type of initial conditions are presented in Sec. IV. We conclude in Sec. V. Appendix A gives some details about the derivation of the transport equations. In Appendix B, we present series solutions of the transport equations that are valid at small $p$. These are used in particular to set appropriate boundary conditions at $p=0$ in the various regimes encountered.

II Transport equations for a quark-gluon system

The analysis, in the framework of kinetic theory, of the evolution of a quark-gluon system towards equilibrium relies on the possibility to describe quark and gluon degrees of freedom in terms of phase space distributions. Color and spin degrees of freedom do not play essential roles in the present discussion and they will be averaged out. We shall denote the color and spin averaged distribution function of gluons with $f(t,{\bf x},{\bf p})$ and that of quarks with $F(t,{\bf x},{\bf p})$ throughout this paper, except in very few cases, such as in eqs. (II.2) or (II.1) below, when a different notation is found more convenient.

In this section we obtain two coupled transport equations that govern the evolution of $f(t,{\bf x},{\bf p})$ and $F(t,{\bf x},{\bf p})$. In a thermal bath of quarks and gluons, the number density of quarks whose masses are much heavier than the temperature $T$ is negligibly small compared to that of light quarks and gluons. We thus only consider the $N_{f}$ flavors of quarks, and their antiparticles, whose masses are smaller than $T$, and take them to be massless for simplicity. Furthermore, we assume that the baryon number density is zero everywhere in the system, and no external forces are exerted on the partons. In this case, quarks and antiquarks have the same distribution due to the $SU(N_{f})$ flavor symmetry and the charge conjugation invariance of QCD. Therefore, one only needs two coupled equations for the quark distribution $F$ and the gluon distribution $f$ to describe the evolution of the system.

Although the number of colors and of flavors are both commonly taken to be $N_{c}=N_{f}=3$ in realistic phenomenological studies of heavy-ion collisions, we will keep here $N_{c}$ and $N_{f}$ as free parameters.

II.1 The Boltzmann equation in the diffusion approximation

The Boltzmann equation

[TABLE]

describes the evolution of the phase space distribution function $f_{\bf p}^{a}$ with the collision term $\mathcal{C}[f_{\bf p}^{a}]$, including all the $2\leftrightarrow 2$ scattering processes in QCD, of the form

[TABLE]

where a short-hand notation $f_{\bf p}^{a}$ is used for the distribution function of different species with the superscript $a$ distinguishing the different particles. Capital letters are used to denote a four-vector, e.g., the four-momentum $P$. Correspondingly, the small and bold letter ${\bf p}$ is used for the three vector, while small ordinary letter $p$ stands for its module. The symbol $\epsilon_{a}$ distinguishes fermions and bosons: $\epsilon_{a}=1$ for bosons and $\epsilon_{a}=-1$ for fermions. In eq. (II.1), the color and spin degrees of freedom of incoming particles $a$ and $b$, and the outgoing particles $c$ and $d$, have been summed over in the squared scattering matrix element $|\mathcal{M}_{cd}^{ab}|^{2}$. The factor $\nu_{a}$ stands for the number of spin $\times$ color degrees of freedom of particle $a$ (which is $2(N_{c}^{2}-1)$ for a gluon and $2N_{c}$ for a quark or an antiquark), and reflects the corresponding averaging of the initial state particle $a$. The factor $s_{cd}$ is a symmetry factor: $s_{cd}$ = $2$ if $c$ and $d$ are identical particles and $s_{cd}=1$ otherwise.

In a pure gluon system, the differential cross-section $gg\leftrightarrow gg$ diverges if the momentum transfer $\bf q$ is much smaller than the momenta of the two scattering gluons. Thus, low momentum transfer or small angle scatterings dominate, which allows us to treat the Boltzmann equation in a diffusion approximation. The kinetic equation then reduces to a Fokker-Planck equation Mueller:1999:Boltzmann ; Blaizot:2013:BEC

[TABLE]

where $\mathcal{J}$ is an effective current that summarizes the effect of the (small angle) collisions. This current is proportional to a logarithmically divergent integral of the form

[TABLE]

where $q_{min}$ is of the order of the screening mass, while $q_{max}$ is of the order of the largest typical momentum in the system (e.g. the temperature if the system is close to equilibrium Arnold:2002:Boltzmann ).

In a quark-gluon system, the small angle scatterings between quarks and gluons are also important. These contribute to two currents: $\mathcal{J}_{g}$ for gluons and $\mathcal{J}_{q}$ for quarks. In addition to the effect of collisions which do not alter the nature of the colliding particles, there are equally important production processes: $q\bar{q}\leftrightarrow gg$, $qg\leftrightarrow qg$ and $\bar{q}g\leftrightarrow\bar{q}g$. These production processes of quarks (gluons) in the scattering of gluons (quarks) with other particles contribute to source terms: $\mathcal{S}_{g}$ for the production of gluons, and $\mathcal{S}_{q}$ for the production of quarks.

By tracking the dominant contributions from all the $2\leftrightarrow 2$ scattering processes between quarks and gluons, as listed in Table 1 of Appendix A, we then obtain two diffusion-like equations

[TABLE]

where the currents are given by

[TABLE]

and the sources by

[TABLE]

with

[TABLE]

Here, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ is the square of the Casimir operator of the color $SU(N_{c})$ group in the fundamental representation. In the following, we shall often refer to the first and second terms on the right hand side of eqs. (II.6) respectively as the diffusion and source terms, although only the contributions proportional to ${\cal I}_{a}$ in the currents correspond truly to diffusion processes.

A few comments on these new equations are in order. First, although we postpone the detailed derivations of eqs. (II.6) to Appendix A, the essential steps and concepts in these derivations can be revealed by focusing on one of the scattering processes, $gq\to gq$, for example. The corresponding diagrams are shown in Fig. 1, and the square of the associated matrix element is

[TABLE]

where the Mandelstam variables are $s=(P+P^{\prime})^{2}$, $t=(P-K)^{2}$ and $u=(P-K^{\prime})^{2}$. There are two types of divergent terms in eq. (II.10) when $q\rightarrow 0$ (small angle approximation). The $u$ channel term ($\sim 1/u$) comes from the square of the second diagram in Fig. 1, while the $t$ channel term ($\sim 1/t^{2}$) comes from the square of the third diagram in Fig. 1. Substituting eq. (II.10) back into eq. (II.1), one finds that the two dominant terms are actually of the same order in the logarithmic approximation †††Here, the difference between the medium-dependent masses of quarks and gluons is neglected, which is valid in the leading logarithmic approximation.. The $t$ channel scattering results in a part of the currents in eqs. (II.7), while the $u$ channel contributes to the sources, eq. (II.8). Repeating the same analysis for all the other $2\leftrightarrow 2$ scattering processes, we obtain eqs. (II.6).

A second comment is that the reduced collision terms in eqs. (II.6) preserve important physical properties of the original kinetic equation, eq. (II.2). For instance, it can be verified that the equilibrium Bose-Einstein distribution for gluons, and Fermi-Dirac distribution for quarks, are still the fixed point solutions to eqs. (II.6) with a temperature given by $T={\cal I}_{a}/{\cal I}_{b}$. Besides, the collision terms in the diffusion form conserve energy, and particle number. We provide an explicit proof for a specified case in the next section.

II.2 The transport equations for spatially homogeneous systems

In the following, we shall study a spatially homogeneous system of quarks and gluons. In this case the spatial dependence of the phase space distribution can be ignored and $D_{t}=\frac{\partial}{\partial t}$. In addition, we assume isotropy of the momentum distributions, which are then solely functions of the modulus of the momentum and of time. We introduce a new time variable

[TABLE]

and denote the derivatives with respect to $\tau$ and $p$ by overdots and primes respectively. Then eqs. (II.6a) and (II.6b) reduce to

[TABLE]

where we have introduced the rescaled currents $J_{g}$ (for gluons) and $J_{q}$ (for quarks), together with the corresponding fluxes $\mathcal{F}_{g}$ and $\mathcal{F}_{q}$:

[TABLE]

and the rescaled source terms

[TABLE]

In the equations above, the integrals $I_{a}$, $I_{b}$ and $I_{c}$ are defined by

[TABLE]

Parton number density $n$, and energy density $\epsilon$, are given in terms of the distribution functions $f$ and $F$ by

[TABLE]

In a similar manner, the entropy density of gluons $s_{g}$ and of quarks $s_{q}$ can be expressed in terms of $f$ and $F$ as

[TABLE]

with the total entropy density of the quark-gluon system given by

[TABLE]

The time evolution of $n$, $\epsilon$ and $s$ can be obtained from eqs. (II.12) and (II.13). The corresponding equations take the following form

[TABLE]

where $s^{+}(p)$ is the non-negative function,

[TABLE]

At this point, it is instructive to discuss the conservation of the total number of partons and of the energy, as well as the increase of the entropy. To do so, one needs to know the behavior of $f$ and $F$ near $p=0$ (the contributions as $p\rightarrow\infty$ to the time derivatives in eqs. (II.22, II.23, II.24) vanish and, therefore, can be dropped). As discussed in Appendix B, two kinds of solutions near $p=0$ are allowed by the transport equations (II.12) and (II.13). For both types of solutions, the boundary terms on the right hand side of eqs. (II.23) and (II.24) always vanish. Therefore, $\dot{\epsilon}=0$ and $\dot{s}\geq 0$. However, $n$ is not conserved with both solutions. For solutions in which $f$ and $F$ are analytic near $p=0$, $n$ is conserved because $\mathcal{F}_{g}$ and $\mathcal{F}_{q}$ vanish at $p=0$. But for solutions of the form

[TABLE]

there is a non-vanishing gluon flux $\mathcal{F}_{g}$ at $p=0$

[TABLE]

This entails a time variation of the number density

[TABLE]

where we have set

[TABLE]

and the coefficient $c_{-1}$ depends only on $\tau$. The non-vanishing gluon flux at $p=0$ reflects the accumulation of gluons of the zero mode, whose number density $N^{0}$ evolves according to

[TABLE]

in order to ensure the overall conservation of the parton number.

In the following, we shall follow Ref. Blaizot:2013:BEC and neglect the mild time dependence of $\mathcal{L}$ in eq. (II.5). In this case eqs. (II.12) and (II.13) are invariant under the following scaling transformation

[TABLE]

with $c>0$. As a result, one can express all momenta in units of $Q_{s}$ (and the same for the chemical potential $\mu$ and the temperature $T$ of Section III) and times in units of $1/Q_{s}$.

III Thermodynamics of QGP with a fixed total parton number

Since only $2\leftrightarrow 2$ processes are included in the collision term of the transport equations, the total parton number is conserved. As a result, in equilibrium, gluons, quarks and antiquarks all have the same chemical potential associated to parton number conservation. The thermal equilibrium distributions are the fixed points of eqs. (II.12) and (II.13), and are of the form

[TABLE]

In the following $T$ and $\mu$ will always refer to as the thermal equilibrium temperature and chemical potential.

The thermodynamic properties of such a QGP are determined by the total energy density and the total parton number density, which are respectively denoted by $\epsilon_{0}$ and $n_{0}$. The under-populated and over-populated systems have very different propertiesBlaizot:2011xf ; Blaizot:2013:BEC . In an under-populated system, the values of $T$ and $\mu<0$ can be obtained by solving the equations

[TABLE]

where $n_{eq}$ and $\epsilon_{eq}$ are obtained by plugging $f_{eq}$ and $F_{eq}$ into eqs. (II.18) and (II.19). In an over-populated system, $n_{0}$ is so large such that no real solution to the above equations exists. The thermal distributions are then given by $f_{eq}$ and $F_{eq}$ with $\mu=0$ and $T$, determined from $\epsilon_{0}$, i.e.,

[TABLE]

The excess gluons form a BEC. The total number of partons with $p>0$, $n_{eq}$, can be calculated from eq. (II.18), and the number density of the condensed gluons is given by

[TABLE]

Let us take for example the system with the total energy and the particle number density

[TABLE]

as obtained from an initial distribution inspired by the color glass picture Blaizot:2011xf (CGC)

[TABLE]

with $f_{0}>0$. The resulting dependence of $\mu$ and $T$ on $f_{0}$ is shown in Fig. 2. The transition from under- to over-population happens at

[TABLE]

Because the production of quarks and antiquarks effectively decreases the number of gluons, larger values of $f_{0c}$ are needed for $N_{f}>0$ than for $N_{f}=0$. For example, $f_{0c}=0.308$ for $N_{c}=3$ and $N_{f}=3$, while $f_{0c}=0.154$ for $N_{c}=3$ and $N_{f}=0$. For $f<f_{0c}$ the system is under-populated. In this case $\mu$ and $T$ can be solved according to eq. (III.34). For $f_{0}>f_{0c}$, the system becomes over-populated. The temperature is then given by eq. (III.35), that is,

[TABLE]

IV Thermalization of the quark-gluon plasma

In this section we study the thermalization of a quark-gluon system whose initial distribution contains only gluons and is of the form given by eq. (III.38). As discussed in the previous section, a BEC is expected to be formed when $f_{0}>f_{0c}$ while when $f_{0}<f_{0c}$ there is no BEC in the equilibrium state. However, we shall show that even in the case $f_{0}<f_{0c}$, when quarks are present, a BEC may appear for a short period of time due to the transient over-population of low momentum gluons. This occurs for $f_{0}>f_{0t}$, where $f_{0t}$ lies in the overlapping region between $\left.f_{0c}\right|_{N_{f}=0}$ and $\left.f_{0c}\right|_{N_{f}>0}$. In the following we study three different patterns of thermalization, each characterized by a specific value of $f_{0}$. In most cases, we take $N_{f}=3$, in which case $f_{0t}\simeq 0.25937<f_{0c}=0.308$.

IV.1 Thermalization with BEC: $f_{0}>f_{0c}$

In Ref. Blaizot:2013:BEC , it is shown that the onset of BEC in a dense system of gluons occurs in a finite time $\tau_{c}$. For the initial condition eq. (III.38), the transition value of $f_{0}$ from under-population to over-population is $\left.f_{0c}\right|_{N_{f}=0}=0.154$, which coincides with the value extracted from eq. (III.39) for $N_{f}=0$. In this subsection, we consider the effects of the quark production on the onset of BEC and manage to follow the evolution of the system, albeit very approximately, beyond $\tau_{c}$. The equilibration process is qualitatively the same for all the over-populated systems with the initial conditions (III.38). We choose $f_{0}=0.4$ and $N_{f}=3$ as a specific example to show the details of how the system evolves into a thermal equilibrium state with BEC.

The formation of BEC starts at a finite time $\tau=\tau_{c}$ when $f$ builds up the $1/p$ tail at small $p$ with the coefficient $c_{-1}$ given by Blaizot:2013:BEC

[TABLE]

As discussed in Appendix B, this is easily understood from the fact that the distribution function at small momentum is accurately described by the classical distribution function,

[TABLE]

with $T^{*}$ and $\mu^{*}_{g}$ time dependent parameters. The onset of BEC corresponds to the vanishing of the effective chemical potential, $\mu^{*}_{g}\to 0$, at which point $f\sim\frac{T^{*}}{p}$. In our numerical simulation, eq. (IV.42) is used to determine the value of $\tau_{c}$ ‡‡‡In our code, $\tau_{c}$ is calculated as the moment when $\left.pf\right|_{p=p_{min}}=T^{*}=\frac{I_{a}}{I_{b}}$ with $p_{min}$ the smallest momentum.. The left panel of Fig. 3 shows how the effective temperature $T^{*}$ keeps decreasing until it eventually approaches the equilibrium temperature $T$. The curve is completely smooth and does not show any indication of the onset of BEC that occurs at $\tau_{c}~{}Q_{s}=0.1708$ (for $f_{0}=0.4$ and $N_{f}=3$). The right panel of Fig. 3 shows the time evolution of the gluon distribution function near $\tau_{c}$, and the approach to the singular behavior, $f(p)\sim 1/p$. The dashed curved are well fitted by the classical distribution (IV.43). Before $\tau_{c}$, $f$ and $F$ are both analytic near $p=0$ and the gluon flux $\mathcal{F}_{g}$ vanishes at $p=0$ (see the left panel of Fig. 4): there is no accumulation of gluons at $p=0$. At $\tau=\tau_{c}$, $f$ becomes singular at $p=0$ but the gluon flux $\left.\mathcal{F}_{g}\right|_{p=0}$ still vanishes according to eqs. (IV.42) and (II.28). As shown in Fig. 5, beyond this moment, the low momentum gluons keep accumulating, as $c_{-1}$ becomes larger than $T^{*}$. Our numerical simulation shows that after $\tau_{c}$ no solutions with vanishing $\left.\mathcal{F}_{g}\right|_{p=0}$ are allowed by the transport equations in (II.12) and (II.13). As discussed in Appendix B, one can find solutions beyond $\tau_{c}$ by providing boundary conditions according to eq. (II.28) (or eq. (B.73)). We have used such solutions to describe the evolution of the system after $\tau_{c}$. Although this procedure ignores important coupling between the condensate and the non-condensate particles, which may alter the details of the dynamics and perhaps the thermalization time scale, it has the advantage of providing a continuous transition to the correct equilibrium state. As shown in the right panel of Fig. 4 and Fig. 5, $\left.\mathcal{F}_{g}\right|_{p=0}$ becomes negative right after $\tau_{c}$, and correspondingly $n$ starts to decreases. This reflects the formation of a BEC, with the number density of condensed gluon, $N^{0}$, increasing according to $\dot{N}^{0}=-\dot{n}$.

The dependence of $\tau_{c}$ on $f_{0}$ can be estimated parametrically at large $f_{0}$. Since $\tau_{c}$ decreases as $f_{0}$ increasesBlaizot:2013:BEC , one needs only study the time evolution of $f$ and $F$ at small $\tau$. This can be done by plugging the linear expansions

[TABLE]

into eqs. (II.12) and (II.13) and keeping terms of $O(\tau^{0})$. We obtain thus

[TABLE]

where

[TABLE]

In the limit $f_{0}\gg 1$ and $p\ll Q_{s}$, we have

[TABLE]

Here, we have dropped the term $\propto\bar{f}^{\prime}_{0}=f_{0}\delta(p-Q_{s})$, which vanishes at $p\ll Q_{s}$. Because the gluon flux vanishes at $p\sim Q_{s}$, $f$ does not change significantly at $p\sim Q_{s}$ and one has $Q_{s}f(Q_{s})\sim Q_{s}~{}f_{0}$. Then $\tau_{c}$ can be estimated as the moment at which $pf$ at small $p$, $f^{3}_{0}Q_{s}^{2}\tau_{c}$, just becomes comparable with $pf$ at $p\sim Q_{s}$, $f_{0}Q_{s}$, which gives

[TABLE]

Since the quark production only contributes a term $\sim-N_{f}\frac{f_{0}^{2}}{p}$ to $\bar{f}_{1}$, eq. (IV.49) is almost independent of $N_{f}$. Thus we do not expect quarks to affect the details of the transition to the BEC when $f_{0}$ is sufficiently large. The parametric behavior (IV.49) is confirmed by the numerical results shown in Fig. 6, and it is actually valid for $f_{0}\gtrsim 1$. Therefore, we conclude that the formation of BEC starts at a time

[TABLE]

for $f_{0}\gtrsim 1$. Note however that for the specific value $f_{0}=0.4$ chosen for the numerical calculations presented in this subsection, $\tau_{c}=0.1708\,Q_{s}^{-1}$, in slight deviation from this relation (which would yield $0.1375\,Q_{s}^{-1}$).

We now consider the effect of quark production on the thermalization process. As we have already mentioned, inelastic processes involving quarks contribute both to the currents and to the source terms in eqs. (II.6). At very early times, the gluon distribution function is approximately given by

[TABLE]

for $p\ll Q_{s}$, and $p$ not too small. The first term on the right hand side of eq. (IV.51) is due to the second part of the current (II.7a), the part proportional to the integral ${\cal I}_{b}$, that drives the increase of the population of soft gluons. The second term is due to the quark production. It acts in the opposite direction, thus hindering the growth of soft gluon modes. However, as shown in Fig. 7, after a short transient period of time, the quark production is peaked at small momenta. This is also confirmed by the plot of the quark flux plotted in Fig. 8: the flux is the largest at small momenta, and continues to increase there all the way till the onset of BEC, and in some cases even beyond the BEC threshold, as revealed by the right hand side of Fig. 8. In this regime, the quark production has no direct effect on the BEC itself. This is because, at small momenta, the outgoing quark current out of a small sphere of radius $p_{0}$ is compensated by the contribution to particle production in that small sphere (i.e. by the source term, as can be verified explicitly by using the small $p$ expansions given in Appendix B, see in particular eq. (• ‣ B) showing that the constant contributions to the current are proportional to ${\cal I}_{c}$ and cancel with the source term, leaving a contribution linear in $p$). This leaves only the gluon current produced by elastic collisions as the source of variation of particle number in the small sphere. And indeed the gluon flux displayed in Fig. 4 is very similar to that obtained for a purely gluonic system (see e.g. Blaizot:2013:BEC ). We have also verified that in the vicinity of the onset, the gluon chemical potential vanishes linearly with ($\tau_{c}-\tau$) within numerical accuracy, as it does in the purely gluonic system.

For the chosen parameters, $N_{f}=3$, $f_{0}=0.4$, the thermal equilibrium can be only achieved after the formation of a gluon BEC. Fig. 9 shows how $f$ and $F$ evolve into thermal distributions after $\tau_{c}$. As we mentioned above, the number of low momentum gluons keeps growing right after $\tau_{c}$. This is a consequence of the relative small rate of quark production (see Fig. 7) and condensate formation $\dot{N}^{0}$, in comparison with the growth rate of low momentum gluons due to the collisions. In the meantime, $\dot{N}^{0}$ increases because of the increase of $c_{-1}$ according to eq. (II.31). The occupation number of low momentum gluons stops growing and starts to decrease at a later time when the condensate formation rate and the quark production rate take over. Note that, as shown in Fig. 7, quark production takes place predominantly at low momentum. The high momentum quark modes are populated by transport. Afterwards, $f$ keeps decreasing while $F$ keeps increasing until the system achieves thermal equilibration. Fig. 10 shows the details about how the number and entropy densities evolve with $\tau$ and eventually reach the predicted values from thermodynamics in Sec. III.

Finally, let us discuss under which conditions the quark production from gluons can be neglected. First, as we have shown in the case $f_{0}\gtrsim 1$, $\tau_{c}$ is (almost) independent of $N_{f}$. On the other hand, for $f_{0}\lesssim 1.0$ quark production delays the onset of BEC and $\tau_{c}$ increases as $N_{f}$ increases. For example, for $f_{0}=0.4$ $\tau_{c}\simeq 0.14~{}Q_{s}^{-1}$ with $N_{f}=0$ and $\tau_{c}\simeq 0.1708~{}Q_{s}^{-1}$ with $N_{f}=3$. This can be easily understood from eq. (IV.46): the production of quarks and antiquarks contributes a negative term $\propto-\frac{N_{f}I_{c}{(0)}}{p}\bar{f}_{0}$ to $\bar{f}_{1}$, which obviously slows down the building-up of the $1/p$ tail of $f$ if $f_{0}$ is not large enough. Second, we observe that the quark production itself slows down the approach to thermalization . To make this statement more quantitative, we define an equilibration time $\tau_{eq}$ by the conditions

[TABLE]

and

[TABLE]

where the values of the above quantities in thermal equilibrium are calculated using $f_{eq}$ and $F_{eq}$ with $T$ and $\mu$ given by eq. (III.41). For $f_{0}=0.4$, we find $\tau_{eq}\simeq 1.1~{}Q_{s}^{-1}$ with $N_{f}=0$ and $\tau_{eq}\simeq 6.4~{}Q_{s}^{-1}$ with $N_{f}=3$. And for $f_{0}=1.0$, we find $\tau_{eq}\simeq 0.86~{}Q_{s}^{-1}$ with $N_{f}=0$ and $\tau_{eq}\simeq 4.8~{}Q_{s}^{-1}$ with $N_{f}=3$. Thus, the presence of quarks increases the thermalization time by typically a factor of 5 (for $N_{f}=3$) (we should keep in mind however that this estimate suffers from the uncertainties related to our very approximate description of the dynamics beyond the onset of BEC).

IV.2 Thermalization without BEC: $f_{0}\leq f_{0t}$

For the initial distribution (III.38), with $f_{0}\leq f_{0t}$, the quark-gluon system will achieve thermal equilibration without the formation of a BEC. Our numerical results verify that the thermal equilibrium temperature $T$ and the negative chemical potential $\mu$ are exactly those predicted by solving eq. (III.34). In those cases, the features of the thermalization process are qualitatively the same for all $f_{0}$. The quarks and antiquarks are produced from the the process $gg\to q\bar{q}$, which causes the gluon number to decrease keeping the total parton number constant. The entropy density of gluons becomes smaller at later times but the total entropy density always increases. Fig. 11 shows the details about how the number and entropy densities of the system with $f_{0}=0.1$ and $N_{f}=3$ evolve into their predicted values in thermal equilibrium. These curves are quite similar to those in Fig. 10, with the noticeable difference that here the parton number $n$ is exactly conserved.

With quark production turned off ($N_{f}=0$), the system of gluons with $f_{0}>\left.f_{0c}\right|_{N_{f}=0}=0.154$ thermalizes with the formation of BEC. As discussed in the previous subsection, the quark production contributes a term $\propto-N_{f}\frac{f_{0}^{2}}{p}\tau$ to $f$ in the early time. For $f_{0t}>f_{0}>\left.f_{0c}\right|_{N_{f}=0}=0.154$, this term is large enough to prevent $f$ from building up the $1/p$ tail near $p=0$, thereby inhibiting the formation of a BEC.

For the same $f_{0}$, the system with $N_{f}\geq 3$ has a lower equilibrium $T$ and a smaller $\mu$ than that with $N_{f}=0$, as shown in Fig. 2. Such a difference causes the under-populated system to thermalize in a different pattern. An example with $f_{0}=0.1$ is shown in Fig. 12. For $N_{f}=0$, the number of low momentum gluons continues to increase until the system achieves thermal equilibrium. For $N_{f}=3$, the occupation number of gluons with $p\lesssim Q_{s}$ first reaches a maximum value which is higher than that in thermal equilibrium. Such an excess of gluons can not be tamed by the quark production until the late stages of equilibration. Let us define two effective chemical potentials

[TABLE]

which are both equal to $\mu$ after the system thermalizes. If $f_{0}<f_{0t}$, $f$ near $p=0$ can be approximated by $f_{eq}$ with $\mu_{g}^{*}/T^{*}<0$. Given $f_{0}$, one can determine the largest value of $\mu_{g}^{*}/T^{*}$ numerically. $f_{0t}$ is defined by the value of $f_{0}$ for which the largest value of $\mu_{g}^{*}/T^{*}$ is zero. At the moment when $\mu_{g}^{*}/T^{*}=0$, $f$ looks like $f_{eq}$ with a vanishing $\mu$ near $p=0$. In the next subsection, we shall show that BEC can be formed due to such a transient excess of gluons in the system with $f_{0}>f_{0t}$ and $N_{f}>0$.

Like in the over-populated case, the quark production delays thermalization. The equilibration time $\tau_{eq}$ is redefined by replacing $n_{g}/n_{geq}$ and $n_{q}/n_{qeq}$ respectively by $\mu^{*}_{g}/\mu$ and $\mu^{*}_{q}/\mu$ in eq. (IV.52). In this definition, the first condition in eq. (IV.52) is a sufficient condition for $f$ and $F$ to be approximately equal to those in thermal equilibrium at small $p$ while the second condition in eq. (IV.53) acts as a constraint to the shape of $f$ and $F$ for the full range of $p$. For $f_{0}=0.1$, we find $\tau_{eq}\simeq 5.5~{}Q_{s}^{-1}$ with $N_{f}=0$ and $\tau_{eq}\simeq 25~{}Q_{s}^{-1}$ with $N_{f}=3$ (again we observe that quark production delays the equilibration time by a factor of $\sim 5$).

IV.3 Thermalization with transient BEC: $f_{0c}>f_{0}>f_{0t}$

A transient BEC can be formed in the under-populated system with $f_{0}>f_{0t}$. Let us choose the system with $f_{0}=0.26$ and $N_{f}=3$ as an example. As shown in the left panel of Fig. 13, $f$ starts to become singular at $p=0$ at $\tau_{c}=0.947~{}Q_{s}^{-1}$. At this moment, the gluon flux $\mathcal{F}_{g}$ still vanishes at $p=0$ because $c_{-1}=T^{*}$ (see eq. (II.29)). However, $c_{-1}$ has a tendency to increase due to the further accumulation of small momentum gluons. Like in the over-populated case, the solution to the transport equations exists after $\tau_{c}$ only if boundary conditions with a non-vanishing $\left.\mathcal{F}_{g}\right|_{p=0}$ are provided. Using the boundary conditions in eq. (B.73) to solve the transport equations, we are able to follow the subsequent evolution of the system. $\mathcal{F}_{g}$ is found to become negative at $p=0$ right after $\tau_{c}$, which is shown in the right panel of Fig. 13. This negative gluon flux reflects the formation of a BEC, and the number density of condensed particles can be calculated from the gluon flux at $p=0$ according to eq. (II.31).

This BEC can only exist for a finite period of time since thermodynamics tell us that the system should evolve into thermal equilibrium without BEC. As shown in the left panel of Fig. 14, $n$ starts to decrease at $\tau=\tau_{c}$, which indicates the formation of BEC. However, $n$ restores its original value after a period $\Delta\tau\simeq 0.35~{}Q_{s}^{-1}$. Afterwards, the solution with vanishing $\left.\mathcal{F}_{g}\right|_{p=0}$ exists again, which describes the subsequent evolution of the system. From that point on $n$ does not change anymore. As expected, $f$, as well as $F$, eventually becomes thermal distributions, which is shown in the right panel of Fig. 14. For an even larger $f_{0}$, the transient BEC exists for a longer time. For example, when $f_{0}=0.28$, we find that BEC starts to form at $\tau_{c}=0.52~{}Q_{s}^{-1}$ and exists for a period of $\Delta\tau\simeq 2.80~{}Q_{s}^{-1}$. In summary, the system with $f_{0c}>f_{0}>f_{0t}$ serves as an example of thermalization with the formation of a transient BEC.

V Discussions

In this paper we have studied the thermalization of a spatially homogeneous quark-gluon plasma, starting from an initial dense system of gluons. Two coupled transport equations for the gluon distribution $f$, and the quark distribution $F$, have been derived using the diffusion approximation of the Boltzmann equation, with the collision term accounting for all possible $2\leftrightarrow 2$ scatterings between quarks and gluons. These transport equations are solved numerically to study how the system evolves from an initial gluon distribution $f_{0}~{}\theta(1-\frac{p}{Q_{s}})$ into a thermalized state of the quark-gluon plasma. We have studied systems with different values of $f_{0}$. $N_{f}$, the number of flavors of quarks that can be taken as massless, is also taken as a free parameter to study the influence of quark production on the formation of BEC and the equilibration process (more precisely, we compare the situation where $N_{f}=3$ to that where $N_{f}=0$). Our main conclusions are

• •

Quark production slows down the growth of $f$ at $p\ll Q_{s}$.

For $N_{f}=0$, a BEC forms for $f_{0}>\left.f_{0c}\right|_{N_{f}=0}=0.154$ in agreement with Ref. Blaizot:2013:BEC . For finite $N_{f}$, there is a range of values of $f_{0}$ larger than $\left.f_{0c}\right|_{N_{f}=0}$ for which quark production hinders the formation of a BEC, and for which the system thermalizes without the formation of a BEC. This occurs for $f_{0c}\leq f_{0}\leq f_{0t}$, where $f_{0t}$ depends on $N_{f}$. We find $f_{0t}\simeq 0.25937$ for $N_{f}=3$.

• •

A transient BEC may develop in intermediate stages prior thermalization.

The critical value $f_{0c}$ characterizing overpopulation depends on $N_{f}$. $f_{0c}=0.308$ for $N_{f}=3$. A BEC is not expected to be formed in equilibrium when $f_{0}<f_{0c}$. However, we find that a transient BEC appears whenever $f_{0c}>f_{0}>f_{0t}$. This is a consequence of the transient excess of low momentum gluons: the growth of low momentum gluon modes is a rapid process, while quark production is relatively much slower. The condensate only exists for a short period of time before quark production eventually takes over and suppresses the excess gluons as the system approaches thermal equilibration.

• •

In the regime of large overpopulation, i.e. for $f_{0}\gtrsim 1$, the formation of BEC occurs at a finite time $t_{c}$ given by the simple formula $t_{c}\sim\frac{1}{(\alpha_{s}f_{0})^{2}}\frac{1}{Q_{s}}$. $t_{c}$ is (almost) independent of $N_{f}$, that is, when $f_{0}$ is large enough, the onset of BEC is not affected by quark production.

• •

Quark production delays thermalization, and quarks are produced dominantly at small momenta. The equilibration time, defined in eqs. (IV.52) and (IV.53), is typically about 5 to 6 times larger for $N_{f}=3$ than that for $N_{f}=0$.

The later observations may have interesting phenomenological consequences, in particular on soft electromagnetic signals Chiu:2012ij , or the elliptic flow Ruggieri:2013bda . However, independently of such potential phenomenological applications, there remain several important theoretical issues that are not addressed in this paper, and that need to be addressed. Like in Refs. Epelbaum:2011pc ; Blaizot:2013:BEC , we only focus on the thermalization of a spatially homogenous non-expanding system. The formation of BEC may also occur in the expanding quark-gluon systemBlaizot:2011xf . It would be of great interest to extend the present work to, say, the boost-invariant $1+1$ dimensional expanding systemMueller:1999:Boltzmann . Moreover, the inelastic processes such as $2\leftrightarrow 3$ are ignored in our transport equations, and it would be important to study how these modify the physical picture that emerges from the present calculation Huang:2013:2to3 . Besides, all the partons are taken as massless and like Ref.s Mueller:1999:Boltzmann ; Venugopalan:2000:Thermalization ; Blaizot:2013:BEC the diffusion approximation is used to simplify the Boltzmann equation. The evolution of the condensates is simply described here by properly added boundary conditions. It would be important to check how reliable those approximations are by a more elaborated investigation on how the low momentum gluons evolve over timecondensate . Finally, the validity of the kinetic description, although widely used in this type of problems, needs to be checked against the statistical classical field simulations, which may be more appropriate at early times Epelbaum:2011pc ; Mueller:Son:2002 . Comparison of the present kinetic approach with the recent studies (see for instance Berges:2013:Turbulent ; Gelis:2013:Isotropization and references therein) would be particularly relevant. We leave all those interesting issues for future studies.

Acknowledgements

We would like to thank F. Gelis for many illuminating discussions. In addition, JPB thanks J. Liao and L. McLerran for collaboration that benefited this work. BW is supported by the Agence Nationale de la Recherche project # 11-BS04-015-01. The research of JPB and LY is supported by the European Research Council under the Advanced Investigator Grant ERC-AD-267258.

Appendix A Diffusion approximation of the collision integral

In this appendix we simplify the collision term of the Boltzmann equation in eq. (II.1) within the diffusion approximationPhysicalKinetics . The squares of the amplitudes for all the $2\leftrightarrow 2$ processes in QCD are listed in Table 1. The momenta of the partons in the final state of these scattering processes are denoted respectively by $K$ and $K^{\prime}$. We only need to keep all the dominant contributions in the limit that the momentum transfer $Q$ is much smaller than the momenta of the two scattering partons, which are denoted respectively by $P$ and $P^{\prime}$. Let us take the $t$ channel dominated processes as an example, in which case $Q=K-P$. In the diffusion limit, the Mandelstam variables reduce to

[TABLE]

with ${\bf v}\equiv{\bf p}/p$ and ${\bf v^{\prime}}\equiv{\bf p^{\prime}}/p^{\prime}$, and

[TABLE]

The corresponding contributions from the $u$ channel scattering can be obtained by simply interchanging $K$ and $K^{\prime}$. The leading contributions to $|{\mathcal{M}}_{cd}^{ab}|^{2}$ in the small angle approximation are given in the third column of Table 1. By plugging the terms proportional to $\frac{s}{t}$ or $\frac{s}{u}$ into the collision term of eq. (II.1), one can easily obtain the source terms

[TABLE]

where we have used the integral

[TABLE]

The terms of $|\mathcal{M}_{cd}^{ab}|^{2}$ proportional to $\frac{s^{2}}{t^{2}}$ and $\frac{s^{2}}{u^{2}}$ in the limit $q\ll p,p^{\prime}$ only contribute to the diffusion terms in the collision term of the transport equations. Let us write

[TABLE]

where the relation of $w^{ab}_{cd}$ to $|{\mathcal{M}}_{cd}^{ab}|^{2}$ can be obtained by referring to eq. (II.1), $\epsilon_{i}=1$ for gluons and $\epsilon_{i}=-1$ for quarks and antiquark. To derive the diffusion terms of the transport equations, one only needs to keep the terms in which the factors in the parentheses $[\cdots]$ on the right hand side of eq. (A.59) vanish in the limit $\bf q\to 0$. In this case, partons $c$ and $d$ can be respectively taken as the same species as $a$ and $b$ §§§Here, we need only to consider the dominant terms from the $t$ channels. There are equal contributions from the $u$ channels if particles $c$ and $d$ are identical particles. However, the sum of the contributions from both channels should be divided by $2$.. Therefore, the diffusion terms describe the diffusion of particle $a$ in the momentum space as a result of scattering off particle $b$. They are different from the source terms, which are proportional to the production rate of particle $b$ of a different species from the scattering parton $a$ with another parton. By expanding the integrand of eq. (A.59) in powers of $q$ and keeping only the first non-vanishing term, we find, after some algebra,

[TABLE]

where the diffusion current for particle $a$ is given by

[TABLE]

with

[TABLE]

Here, $|{\mathcal{M}}_{ab}^{ab}|^{2}_{diff}$ are the terms proportional to $\frac{s^{2}}{t^{2}}$ in the third column of Table 1. To simplify $\mathcal{J}^{a}$, we need to evaluate

[TABLE]

and

[TABLE]

Here, we have assumed that $f^{b}_{\bf p}=f^{b}_{-\bf p}$ in order to get $J^{ab}$ in the last line in the above equation. $\mathcal{J}^{a}$ in eq. (II.7) is obtained by summing $J^{ab}$ over $b$ with the coefficient given by that of the corresponding term proportional to $\frac{s^{2}}{t^{2}}$ in the third column of Table 1.

Appendix B Series solutions and boundary conditions to the transport equations

As discussed in the main text, there are two types of solutions of the transport equations, characterized by the behavior of the gluon distribution near the origin $p=0$: either $f(p=0)$ is a finite constant, or $f(p\rightarrow 0)\sim 1/p$. In order to analyze further these solutions, we set

[TABLE]

where the coefficients $c_{n}$ and $d_{n}$ can be determined from the transport equations in eqs. (II.12) and (II.13) with $I_{a}$, $I_{b}$ and $I_{c}$ functions of $\tau$. One then finds that there are only two types of solutions allowed by the transport equations:

• •

$f$ is analytic at $p=0$.

In this case, we have

[TABLE]

and

[TABLE]

with $c_{0}=f(\tau,0)$ and $d_{0}=F(\tau,0)$.

In the limit $c_{0}\gg 1$, the radius of convergence of the above series solution shrinks to zero. In this case, we find

[TABLE]

with $T^{*}\equiv\frac{I_{a}}{I_{b}}$. The leading terms in $c_{0}$ at each order in $p$ can be resummed and, thus, we obtain

[TABLE]

This is the classical distribution function, with an effective chemical potential given $\mu^{*}_{g}\equiv-T^{*}/c_{0}$. The above resummed solution is very useful for understanding the evolution of the quark-gluon system close to $\tau_{c}$Blaizot:2013:BEC .

• •

$f$ is singular at $p=0$.

In this case, we get

[TABLE]

and

[TABLE]

To solve the transport equations in eqs. (II.12) and (II.13), one needs two initial conditions and four boundary conditions. In our code, we use the following boundary conditions

[TABLE]

The explicit Euler method is used for time integration and only $c_{-1}$ at the current time step is needed for the calculation of $f$ and $F$ at the next time step.