Effective kinetic description of the expanding overoccupied Glasma

TL;DR

This study numerically investigates the kinetic evolution of an overoccupied, expanding Glasma, revealing gluon and quark distribution behaviors that resemble classical simulations and suggest universal self-similar scaling.

Contribution

It demonstrates that a kinetic Boltzmann approach captures key features of Glasma evolution, including gluon cascades and self-similar quark distributions, aligning with classical-statistical results.

Findings

Gluon distributions show infrared enhancement similar to Bose-Einstein condensation.

Quark distributions exhibit self-similar scaling behavior.

Results challenge traditional thermalization scenarios.

Abstract

We report on a numerical study of the Boltzmann equation including $2\leftrightarrow 2$ scatterings of gluons and quarks in an overoccupied Glasma undergoing longitudinal expansion. We find that when a cascade of gluon number to the infrared occurs, corresponding to an infrared enhancement analogous to a transient Bose-Einstein condensate, gluon distributions qualitatively reproduce the results of classical-statistical simulations for the expanding Glasma. These include key features of the distributions that are not anticipated in the "bottom-up" thermalization scenario. We also find that quark distributions, like those of gluons, satisfy self-similar scaling distributions in the overoccupied Glasma. We discuss the implications of these results for a deeper understanding of the self-similarity and universality of parton distributions in the Glasma.

Equations

$$f(\tau ,\mathbf{p})=\frac{(2\pi {)}^3}{2(N_c^2-1)}\frac{d{N}_{\text{gluon}}}{d^{3}x{d}^{3}p}.$$

$$F(\tau ,\mathbf{p})=\frac{(2\pi {)}^3}{2N_c{N}_f}\frac{d{N}_{\text{quark}}}{d^{3}x{d}^{3}p}.$$

$$\left(\frac{\partial }{\partial \tau }-\frac{p_z}{\tau }\frac{\partial }{\partial p_z}\right)f(\tau ,\mathbf{p})$$

$$\left(\frac{\partial }{\partial \tau }-\frac{p_z}{\tau }\frac{\partial }{\partial p_z}\right)F(\tau ,\mathbf{p})$$

$$C_{\text{gluon}}[f,F]$$

$$C_{\text{quark}}[F,f]$$

$${\mathbf{J}}_{\text{g}}$$

$${\mathbf{J}}_{\text{q}}$$

$$S_{\text{g}}$$

$$S_{\text{q}}$$

$$I_a(\tau )=\int \frac{d^{3}p}{(2\pi {)}^3}\left[N_{c}f(1+f)+N_{f}F(1-F)\right],$$

$$I_b(\tau )=2\int \frac{d^{3}p}{(2\pi {)}^3}\frac{1}{p}\left[N_{c}f+N_{f}F\right],$$

$$I_c(\tau )=\int \frac{d^{3}p}{(2\pi {)}^3}\frac{1}{p}\left[f+F\right].$$

$$\mathcal{L}={\int }_{q_{\text{min}}}^{q_{\text{max}}}\frac{dq}{q}=\log \frac{q_{\text{max}}}{q_{\text{min}}}.$$

$$m_D=\sqrt{2g^2{I}_b(\tau )}$$

$$⟨p_{\perp }^2⟩=\frac{\int d^{3}p{p}_{\perp }^{2}f(\tau ,p_{\perp },p_z)}{\int d^{3}pf(\tau ,p_{\perp },p_z)},$$

$$\tilde{f}(\tau ,\mathbf{p})=\frac{f(\tau ,\mathbf{p})}{f_0},$$

$$\tilde{\tau }=\frac{1}{4\pi }(N_c{g}^2{f}_0{)}^2{\mathcal{L}}_0\tau ,$$

$$\left(\frac{\partial }{\partial \tilde{\tau }}-\frac{p_z}{\tilde{\tau }}\frac{\partial }{\partial p_z}\right)\tilde{f}(\tilde{\tau },\mathbf{p})=\tilde{\mathcal{L}}{\nabla }_{\mathbf{p}}\cdot \left[{\tilde{I}}_a{\nabla }_{\mathbf{p}}\tilde{f}(\tilde{\tau },\mathbf{p})+\frac{\mathbf{p}}{p}{\tilde{I}}_b\tilde{f}(\tilde{\tau },\mathbf{p})\left(f_0^{-1}+\tilde{f}(\tilde{\tau },\mathbf{p})\right)\right]+\frac{C_F{N}_f}{N_c^2{f}_0}\tilde{S},$$

$$\left(\frac{\partial }{\partial \tilde{\tau }}-\frac{p_z}{\tilde{\tau }}\frac{\partial }{\partial p_z}\right)F(\tilde{\tau },\mathbf{p})=\frac{C_F}{N_c}\tilde{\mathcal{L}}{\nabla }_{\mathbf{p}}\cdot \left[{\tilde{I}}_a{\nabla }_{\mathbf{p}}F(\tilde{\tau },\mathbf{p})+\frac{\mathbf{p}}{p}{\tilde{I}}_{b}F(\tilde{\tau },\mathbf{p})\left(1-F(\tilde{\tau },\mathbf{p})\right)\right]-\frac{C_F^2}{N_c^2}\tilde{S},$$

$$\tilde{S}=\tilde{\mathcal{L}}{\tilde{I}}_c\frac{1}{p}\left[F(f_0^{-1}+\tilde{f})-\tilde{f}(1-F)\right].$$

$${\tilde{I}}_a(\tilde{\tau })=\int \frac{d^{3}p}{(2\pi {)}^3}\left[\tilde{f}(f_0^{-1}+\tilde{f})+\frac{N_f}{N_c{f}_0^2}F(1-F)\right],$$

$${\tilde{I}}_b(\tilde{\tau })=2\int \frac{d^{3}p}{(2\pi {)}^3}\frac{1}{p}\left[\tilde{f}+\frac{N_f}{N_c{f}_0}F\right],$$

$${\tilde{I}}_c(\tilde{\tau })=\int \frac{d^{3}p}{(2\pi {)}^3}\frac{1}{p}\left[\tilde{f}+\frac{1}{f_0}F\right],$$

$$\tilde{\mathcal{L}}(\tilde{\tau })=\frac{\mathcal{L}(\tau )}{{\mathcal{L}}_0}.$$

$$m_D^2=2N_c{g}^2{f}_0{\tilde{I}}_b(\tilde{\tau }).$$

$$n(\tau )=n_{\text{g}}(\tau )+n_{\text{q}}(\tau ),$$

$$\mathcal{E}(\tau )={\mathcal{E}}_{\text{g}}(\tau )+{\mathcal{E}}_{\text{q}}(\tau ),$$

$$n_{\text{g}}(\tau )=2(N_c^2-1)\int \frac{d^{3}p}{(2\pi {)}^3}f(\tau ,\mathbf{p}),$$

$$n_{\text{q}}(\tau )=4N_c{N}_f\int \frac{d^{3}p}{(2\pi {)}^3}F(\tau ,\mathbf{p}),$$

$${\mathcal{E}}_{\text{g}}(\tau )=2(N_c^2-1)\int \frac{d^{3}p}{(2\pi {)}^3}pf(\tau ,\mathbf{p}),$$

$${\mathcal{E}}_{\text{q}}(\tau )=4N_c{N}_f\int \frac{d^{3}p}{(2\pi {)}^3}pF(\tau ,\mathbf{p}).$$

Full text

I Introduction

There are several remarkable and puzzling features of the early time nonequilibrium evolution of quark-gluon matter in ultrarelativistic heavy-ion collisions. A systematic study of its properties is feasible in the limit where the QCD coupling is very weak and the occupancies of gluons are correspondingly large such that their product $g^{2}\,f_{g}\sim\mathcal{O}(1)$. In this limit, the spacetime evolution of the overoccupied matter–the Glasma Kovner:1995ja ; Lappi:2006fp –can be studied employing real-time classical-statistical simulations Krasnitz:1998ns ; Lappi:2003bi of the Yang-Mills equations. Recent large scale 3+1-dimensional Yang-Mills simulations of the expanding Glasma clearly demonstrate that it flows to a nonthermal turbulent fixed point characterized by a self-similar gluon distribution Berges:2013eia ; Berges:2013fga . Rather unexpectedly, the simulations showed that key properties of the gluon distribution were identical to those anticipated in the “bottom-up” thermalization scenario of Baier et al. Baier:2000sb .

The results of the classical-statistical simulations are unexpected because gluon dynamics in the expanding Glasma shows no trace of late time plasma instabilities–these should be present in weak coupling kinetic frameworks Mrowczynski:1993qm ; Kurkela:2011ti ; Mrowczynski:2016etf . Further, the ratio of the longitudinal pressure to the transverse pressure $P_{L}/P_{T}$ deviates from the expectation from kinetic theory–including as well, that from the bottom-up scenario. Finally, and most remarkably, the properties of the nonthermal Glasma fixed point are, in a wide inertial range of momenta, identical to those of self-interacting scalar fields prepared with the same geometry Berges:2014bba . Put together, these observations pose a challenge for straightforward kinetic descriptions of the expanding Glasma Berges:2015ixa .

In this paper, we will investigate the nonequilibrium early time evolution of quark and gluon matter in the longitudinally expanding geometry of an ultrarelativistic heavy-ion collision by numerically solving the Boltzmann equation with an overoccupied Glasma initial condition. Such a study was previously performed111There are studies of the Boltzmann equation with an expanding geometry but they have not included yet the $(1+f)$ Bose enhancement factor that is essential to match the classical Yang-Mills simulations to kinetic theory Ruggieri:2015tsa ; Greif:2016jeb . by Kurkela and Zhu Kurkela:2015qoa , who implemented the effective kinetic theory of Arnold, Moore and Yaffe Arnold:2002zm including both elastic 2$\leftrightarrow$2 scattering and soft splitting processes. They however restricted their simulations to include only gluons. Further, their focus was to follow the evolution of the Glasma all the way to thermalization a la the bottom-up scenario.

Our interest here is primarily the classical stage of the Glasma and we will perform simulations for a wider range of couplings than those presented in Kurkela:2015qoa . Firstly, by including only elastic 2$\leftrightarrow$2 processes, we will investigate how well these kinetic studies reproduce the self-similar scaling behavior of the gluon distribution, that has been obtained by the classical-statistical simulations. We will employ the small angle approximation to the Boltzmann equation whereby, as shown by Landau landau1936kinetic ; lifshitz1981physical , the collision integral is approximated by a Fokker-Planck like diffusion term222The Landau kinetic equation was first applied to the Glasma by Mueller Mueller:1999pi and simulated numerically in Bjoraker:2000cf . However in Refs. Mueller:1999pi ; Bjoraker:2000cf , the Bose enhancement factor, important for the description of overoccupied gluons, was not taken into account.. We will demonstrate that a flow of gluon number to the infrared (IR) generates an IR enhancement that often interpreted as a transient Bose-Einstein condensate (BEC) in the literature. We will present results for both the longitudinal and transverse momentum distributions as well as the behavior of $P_{L}/P_{T}$, demonstrating the influence of this IR enhancement. We should emphasize at the outset that the IR enhancement is unlikely to capture key features of BEC formation as studied and quantified in detail for identical overoccupied conditions in scalar theories off-equilibrium Berges:2014bba ; Berges:2015ixa .

The generation of such an IR enhancement is unsurprising for a system of overoccupied gluons, as argued in a number of papers Blaizot:2011xf ; Blaizot:2013lga ; Xu:2014ega . What is surprising is that, when the BEC becomes manifest, our results reproduce important qualitative features of the classical-statistical simulations of the Glasma. If inelastic scattering processes played a large role at these early times, a BEC would never form and would not influence the dynamics of the Glasma. Our results therefore suggest that, at least for the range of weak couplings studied, inelastic processes may not be sufficiently strong at early times to hinder the formation of a transient BEC.

We emphasize that our framework does not provide a first principles solution to some of the puzzling dynamics of the Glasma that we listed. It does not contain either $2\leftrightarrow 3$ scattering or the plasma instabilities that should be there in a purely kinetic picture. Further, it is unclear how one defines a BEC of gluons in QCD. Our result should instead be thought of as an effective description of the classical-statistical numerical simulations for the range of couplings in those studies. Both inelastic scattering processes and late time plasma instabilities, whose dynamics is in principle contained in the classical-statistical simulations, appear to be suppressed in the early-time evolution of the overoccupied Glasma. Our results therefore argue for a deeper understanding and characterization of the IR dynamics in the Glasma. Key features of this dynamics are captured by the spatial string tension, which has a characteristic time dependence Mace:2016svc that is universal Berges:2017igc . It has been observed that the spatial string tension sets the scale for off-equilibrium sphaleron transitions Mace:2016svc ; this latter quantity is essential for dynamical simulations Mueller:2016ven ; Mace:2016shq of the Chiral Magnetic effect.

The other novel feature of our study is the kinetic description of quark distributions in the expanding Glasma. The Landau kinetic equations were studied for an isotropic distribution of overoccupied gluons along with quarks and anti-quarks Blaizot:2014jna ; our work extends this study to the expanding Glasma. We will show that quark distributions in the Glasma also obey self-similar scaling distributions with the same dynamical scaling exponents as those satisfied by the gluons. In an accompanying paper, we applied the insights from this study to investigate photon production in the Glasma relative to the quark-gluon plasma Berges:2017eom . For completeness, we will study chemical equilibration in the Glasma–while interesting, the results here for the relative population of quarks and gluons will be modified by the collinear gluon splitting processes that also influence kinetic equilibrium. Their additional contribution will be explored in subsequent work.

The paper is organized as follows. In the next section, we will outline the kinetic equations for the coupled system of gluons, quarks and anti-quarks in the Glasma. Our discussion here is a straightforward extension of that in Blaizot:2014jna . Numerical results are discussed in section 3 for the Glasma with no quarks or anti-quarks ($N_{f}=0$). We consider the time evolution of the longitudinal and transverse hard scales, the ratio of the longitudinal and transverse pressure, and the self-similarity of the gluon distribution. For very weak couplings, qualitative agreement is found with the classical-statistical simulations, while significant deviations are observed for larger couplings. In section 4, we extend our study to finite $N_{f}$. The time evolution of the quark distributions is studied for a range of initial occupancies. As noted, the self-similar behavior that characterizes the gluon distribution holds for quark distributions as well. The approach to chemical equilibration is studied and is shown to be sensitive to the initial quark occupancy. The final section contains a brief summary of our results and a discussion of future work. The numerical procedure is outlined in Appendix A and the features of the BEC formed are discussed in Appendix B.

II Kinetic equation for 2$\leftrightarrow$2 elastic scatterings in the Glasma

II.1 Small angle approximation for the Boltzmann equation

We begin by defining the gluon distribution function as a number density in phase space averaged over all polarization and color states:

[TABLE]

Similarly, for $N_{f}$-flavor massless and unpolarized quarks, the distribution function is defined as

[TABLE]

We will assume that the system is boost-invariant in the longitudinal direction and uniform in the transverse direction, so that $f$ and $F$ are functions of $\tau=\sqrt{t^{2}-z^{2}}$ and ${\boldsymbol{p}}$ alone. Furthermore, we assume that anti-quarks have the same distribution as quarks. Each of these distributions satisfies a Boltzmann equation, defined as

[TABLE]

where $-\frac{p_{z}}{\tau}\frac{\partial f}{\partial p_{z}}$ is the drift term describing the effect of the longitudinal expansion. For the collision terms $C_{\text{gluon}}[f,F]$ and $C_{\text{quark}}[F,f]$, we consider only 2$\leftrightarrow$2 elastic scatterings. Since gluons, quarks and anti-quarks scatter off each other, the collision terms of each of the equations depend on the distributions of the other, which makes Eqs. (3) and (4) a coupled set of equations. In the small angle approximation, the collision terms can be expressed as the sum of a diffusion term and a source term Blaizot:2014jna ,

[TABLE]

Here ${\boldsymbol{J}}_{\text{g}}$ and ${\boldsymbol{J}}_{\text{q}}$ denote particle currents in momentum space and are given by

[TABLE]

The source terms $S_{\text{g}}$ and $S_{\text{q}}$ are respectively

[TABLE]

where $C_{F}=(N_{c}^{2}-1)/(2N_{c})$ is the quadratic Casimir in the fundamental representation of SU$(N_{c})$. In obtaining these expressions, it has been assumed that $f(\tau,{\boldsymbol{p}})=f(\tau,-{\boldsymbol{p}})$ and $F(\tau,{\boldsymbol{p}})=F(\tau,-{\boldsymbol{p}})$. We have further defined the integrals

[TABLE]

[TABLE]

and

[TABLE]

The integrand of $I_{a}$ is proportional to the density of possible scatterers that are enhanced or suppressed by the Bose or the Pauli factor respectively. The integral $I_{b}$ is related to the Debye mass scale $m_{D}$ as $m_{D}^{2}=2g^{2}I_{b}$. The Coulomb logarithm $\mathcal{L}$ is a divergent integral that is regularized by cutoffs as

[TABLE]

This IR divergence originates from the long range nature of the interaction and is regularized by the medium mass. We will use the Debye mass scale

[TABLE]

as the IR cutoff. Since $m_{D}$ is time dependent, so too is the Coulomb logarithm.

The ultraviolet (UV) divergence, in contrast, is not physically meaningful since it arises from the small momentum expansion performed inside the momentum integrals that give the above form of the kinetic equations. The original collision integral does not possess this UV divergence because the distribution functions have finite support in momentum space. It is therefore natural to use the typical momentum scale of the distribution functions as the UV cutoff. Since the typical transverse momentum scale is larger than the longitudinal one in a longitudinally expanding Glasma, we will use the square root of the mean $p_{\!\perp}^{2}$,

[TABLE]

as the UV cutoff for the Coulomb logarithm. In the overoccupied weak coupling regime of interest in this study, this scale is nearly constant and close to the saturation scale $Q_{s}$. We will therefore simply choose the UV cutoff to be $Q_{s}$, except in Sec. IV.2.

II.2 Rescaling for an overoccupied Glasma

In the following, we will exploit the fact that the Glasma is highly overoccupied at early times. To simplify further analysis, we will rescale the gluon distribution as

[TABLE]

with $f_{0}$ being the maximum value of the initial distribution $f(\tau_{0},{\boldsymbol{p}})$. We further introduce the rescaled time $\tilde{\tau}$,

[TABLE]

with $\mathcal{L}_{0}$ denoting the initial value of the Coulomb logarithm, $\mathcal{L}_{0}=\mathcal{L}(\tau_{0})$, determined by the initial distributions. Since the quark occupancy is always smaller than one, we will not rescale it.

In terms of the rescaled quantities, the kinetic equations can be expressed as

[TABLE]

for gluons, and

[TABLE]

for quarks, with the rescaled source term

[TABLE]

The corresponding rescaled integrals are

[TABLE]

[TABLE]

[TABLE]

and we have introduced a normalized Coulomb logarithm defined as

[TABLE]

In terms of the rescaled integral $\tilde{I}_{b}$, the Debye mass is expressed as

[TABLE]

Hereafter, we shall fix $N_{c}=N_{f}=3$. The relevant physical parameters for this system are the coupling $g$, the initial occupancy $f_{0}$, and the initial time $\tilde{\tau}_{0}$. To ensure that a kinetic approach based on a perturbative collision kernel is valid, both $g$ and $g^{2}f_{0}$ must be sufficiently small.

Comparing, Eqs. (19) and (20), we can make the following observations.

• •

Since $C_{F}/N_{c}=4/9<1$, the effect of diffusion is weaker for quarks than gluons.

• •

If the initial state is highly occupied by gluons, i.e. $f_{0}\gg 1$, the source term is much smaller than the diffusion term for gluons, while for quarks the source term is always as important as the diffusion term at least parametrically.

II.3 Conservation laws

One can easily confirm that the collision terms, as approximated by Eqs. (5) and (6), conserve both the total particle number density

[TABLE]

and the total energy density

[TABLE]

where the particle number density and the energy density, respectively, of each species are

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In the presence of the expansion term, the conservation laws are modified to

[TABLE]

and

[TABLE]

where $\mathcal{P}_{L}$ is the longitudinal pressure given by the sum of

[TABLE]

and

[TABLE]

For later convenience, we also introduce the transverse pressure

[TABLE]

II.4 Initial conditions

For a given initial distribution $f(\tau_{0},{\boldsymbol{p}})$, the initial time $\tau_{0}$ which is consistent with that initial distribution can be specified by the following argument. The number density of hard gluons produced immediately after a heavy-ion collision can be expressed as Mueller:1999fp

[TABLE]

where $c$ is the gluon liberation coefficient, which can be estimated by using classical Yang-Mills simulations of the Glasma Krasnitz:2000gz ; Krasnitz:2003jw ; Lappi:2003bi . We employ the value $c=1.1$ given in Lappi:2007ku . As long as elastic scattering dominates amongst hard gluons, the number density decreases in time as $1/\tau$ due to the longitudinal expansion. Since the gluon number density is dominated by hard gluons at early times, the expression (38) should agree with $n_{g}$ given by (29). By this condition, we can find the initial time $\tau_{0}$ that is consistent with the initial distribution $f_{0}({\boldsymbol{p}})=f(\tau_{0},{\boldsymbol{p}})$ as

[TABLE]

The rescaled time that corresponds to this initial time is

[TABLE]

III Numerical results I. $N_{f}=0$

In this section, we will show numerical results for a pure gluon system ($N_{f}=0$). We will employ a family of initial gluon distributions,

[TABLE]

where $p_{\!\perp}=\sqrt{p_{x}^{2}+p_{y}^{2}}$ is the transverse momentum, the parameter $\xi_{0}$ characterizes the anisotropy of the initial distribution and $n_{0}=g^{2}f_{0}$ specifies the overoccupancy. For the kinetic description to be valid, $n_{0}$ must be sufficiently small, as it controls the effective interaction strength.

Substituting this initial distribution into Eq. (40), one obtains the initial time

[TABLE]

where the initial value of the Coulomb logarithm is

[TABLE]

In the following, we will mostly use $n_{0}=0.1$. The corresponding numerical values of the initial time are $Q_{s}\tilde{\tau}_{0}\simeq 0.33$, 0.74 and 1.7 for $\xi_{0}=1$, 2 and 4, respectively. For simplicity, we shall typically use the initial time of $Q_{s}\tilde{\tau}_{0}=1$ except in Sec. IV.2.

To solve the kinetic equation in Eq. (19) numerically, we use $p=\sqrt{p_{\!\perp}^{2}+p_{z}^{2}}$ and $\kappa=\cos\theta=p_{z}/p$ as the momentum variables. For $p$, we discretize the range $p_{\text{min}}\leq p\leq p_{\text{max}}$ nonuniformly into $N_{p}$ grid points. For the angular variable $\kappa$, it is sufficient to consider the range $0\leq\kappa\leq 1$, which is divided into $N_{\kappa}$ grid points. Unless otherwise noted, we use $p_{\text{min}}/Q_{s}=10^{-2}$, $p_{\text{max}}/Q_{s}=8$, $N_{p}=500$ and $N_{\kappa}=256$. We have confirmed that our results are insensitive to changes of these parameters. The details of the numerical method are discussed in Appendix A.

III.1 Self-similar evolution

We will begin with a brief review of the classical regime (of gluon occupancies greater than unity) in the bottom-up thermalization scenario ($Q_{s}^{-1}\ll\tau\ll\alpha_{s}^{-3/2}Q_{s}^{-1}$) Baier:2000sb . In this time range, the number of gluons is dominated by hard gluons whose transverse momenta are the order of the saturation scale $Q_{s}$, and quarks and gluons scatter at small angles with typical momentum transfers $\sim m_{D}$. The initial density of hard gluons is parametrically $Q_{s}^{3}/\alpha_{s}$, and it decreases in time due to the longitudinal expansion as

[TABLE]

The numerical coefficient for this relation is given by Eq. (38); however, it is not necessary for this qualitative discussion. The Debye mass can be estimated as

[TABLE]

Since the exchanged momentum is much smaller than the saturation scale, the typical transverse momentum of hard gluons stays around $Q_{s}$. The small angle scatterings primarily modify the typical longitudinal momentum, which we call $p_{z}$, by providing random kicks of order $m_{D}$. Since these random and incoherent scatterings can be regarded as a diffusion process, the typical longitudinal momentum obeys

[TABLE]

Here, $N_{\text{col}}$ is the number of the collisions. Using $f_{\text{hard}}\sim n_{\text{hard}}/(Q_{s}^{2}p_{z})$, one can estimate it as

[TABLE]

where $(1+f_{\text{hard}})$ is the Bose enhancement factor and $\sigma\sim\alpha_{s}^{2}/m_{D}^{2}$ is the cross section. By combining (44), (46) and (47), we obtain

[TABLE]

and

[TABLE]

These results can be summarized by the scaling relation

[TABLE]

where $f_{S}$ is a scaling function. As noted previously, this scaling law has been confirmed by large scale classical-statistical real-time lattice simulations Berges:2013eia ; Berges:2013fga . The keys to the above derivation are the longitudinal expansion, the use the diffusion relation Eq. (46) and the Bose enhancement factor in Eq. (47). All of these features are included in the kinetic equation with the small angle approximation Eq. (19). Indeed, the diffusion constant $I_{a}$ can be interpreted as $dN_{\text{col}}/d\tau$. Therefore we anticipate that bottom-up thermalization in the overoccupied Glasma is well described by the kinetic equation (19).

In the above discussion, the soft splitting processes are not considered. Although those number-changing processes are naively higher order in the coupling, they give the leading order contributions because of the collinear enhancement Arnold:2002zm . However, their main effect is producing soft gluons that have transverse momentum much smaller than $Q_{s}$, and the self-similar scaling behavior shown by the hard gluons is not affected by them. In the bottom-up thermalization scenario, the soft splitting processes play a crucial role in the later stage $\tau\gg\alpha_{s}^{-3/2}Q_{s}^{-1}$ Baier:2000sb .

In Fig. 1, the Debye mass squared is plotted as a function of the rescaled time. Results for several values of the coupling, $g=10^{-3}$, $10^{-2}$ and $10^{-1}$, are compared. Other parameters are $n_{0}=0.1$, $\xi_{0}=2$ and the initial time333These parameters correspond to the initial time $Q_{s}\tau_{0}\simeq 52$ in the original time variable. $Q_{s}\tilde{\tau}_{0}=1$. Since we fix $n_{0}$, the Debye mass is independent of the coupling in the classical scaling regime. The $1/\tau$ decrease in Eq. (45) can be confirmed clearly in the simulations.

The time evolution of the rescaled distribution $\tilde{f}=f/f_{0}$ evaluated at $p_{z}=0$ is plotted as a function of the transverse momentum in the left panel of Fig. 2. The parameters used for this computation are $g=10^{-3}$, $n_{0}=0.1$ ($f_{0}=10^{5}$), $\xi_{0}=2$ and $Q_{s}\tilde{\tau}_{0}=1$. We note that in the infrared, the distribution shows $1/p_{\!\perp}$ behavior. This result reproduces qualitatively444In detail, we note that the shape of the distribution near the hard scale is very similar to that observed for the scalar theory discussed at length in Berges:2015ixa . In the classical-statistical simulations of the Glasma, this “bump” was not seen for the times studied. We thank Juergen Berges for this observation. those for $p_{\!\perp}$ distributions from the classical-statistical simulations of the expanding Glasma Berges:2013eia ; Berges:2013fga . In this kinetic framework, the emergence of the $1/p_{\!\perp}$ behavior in the Glasma is a consequence of the onset of the transient Bose-Einstein condensate Blaizot:2013lga ; Blaizot:2014jna . The onset of the BEC within our approximation is discussed in Appendix B. For a proper description of the IR soft region, number-changing inelastic processes are essential Blaizot:2016iir ; further, it likely requires that one goes beyond a conventional kinetic framework and considers resummations along the lines discussed in Orioli:2015dxa . Our primary interest here lies in the universal scaling behavior that the hard modes show.

In the right panel of Fig. 2, we plot the transverse momentum distribution multiplied by the scaling factor $(Q_{s}\tilde{\tau})^{2/3}$. For hard momenta p_{\!\perp}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}Q_{s}, the curves at different times nicely overlap confirming the scaling behavior suggested by Eq. (50). Figure 3 shows similar plots for the second moment of the transverse distribution, $p_{\!\perp}^{2}\tilde{f}(\tilde{\tau},p_{\!\perp},0)$. From these plots, one can confirm that the typical transverse momentum scale is of the order of $Q_{s}$. In Fig. 4, the time evolution of the longitudinal momentum distribution evaluated at $p_{\!\perp}=Q_{s}$ is plotted. The left panel shows the original distribution, while the right panel shows the distribution multiplied by the scaling factor $(Q_{s}\tilde{\tau})^{2/3}$ as a function of the rescaled momentum $(Q_{s}\tilde{\tau})^{1/3}p_{z}$. Again, the rescaled distributions for vastly different times lie on top of each other. These results confirm that the nonequilibrium evolution of weakly coupled, overoccupied gluons in a longitudinally expanding system follows the bottom-up thermalization scenario.

We also numerically compute the perturbative expressions for the hard scale observables Berges:2013fga ,

[TABLE]

and

[TABLE]

In the scaling regime, the longitudinal hard scale is expected to behave as $\Lambda_{L}^{2}(\tau)\sim(Q_{s}\tau)^{-2/3}$, while the transverse hard scale stays nearly constant. Their respective behavior is confirmed by the numerical results shown in Fig. 5, where the hard scales are plotted as a function of time for different values of the coupling. For the longitudinal hard scale, the free streaming behavior, $\Lambda_{L}^{2}(\tau)\sim(Q_{s}\tau)^{-2}$, is seen at early times. At later times, the system is driven to the classical scaling regime indicated by the temporal exponent $-2/3$ for sufficiently weak coupling, g\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}10^{-2}. On the other hand, results for the larger couplings (smaller occupancies) clearly deviate from the anticipated classical scaling behavior.

The ratio of the longitudinal pressure to the transverse pressure is a measure of the bulk anisotropy of the system. As shown in Fig. 6, the ratio decreases in time because of the longitudinal expansion of the system. For larger couplings, it turns to increase at later times and approaches unity. For weak couplings, the system falls into the scaling regime. If the pressure is dominated by hard modes, this ratio should show the temporal power law behavior $\sim(Q_{s}\tau)^{-2/3}$ in the scaling regime. However, as seen in Fig. 6, the curves clearly deviate from the scaling behavior $(Q_{s}\tau)^{-2/3}$ – more significantly than that observed for the longitudinal hard scale. The same tendency was noted in the classical-statistical simulations; this is because the longitudinal pressure is not dominated by the hard modes and soft modes contribute significantly Berges:2015ixa .

To check the robustness of the classical scaling behavior, we compare different initial distributions for the longitudinal hard scale in Fig. 7. For the initial distributions given by Eq. (41), we compute results for anisotropy parameter values of $\xi_{0}=1$, 2 and 4. Furthermore, results for a different functional form of the initial distribution,

[TABLE]

denoted as Box IC, are compared to the results from distributions with the initial conditions in Eq. (41). The other parameters in the computation are set to $g=10^{-2}$, $n_{0}=0.1$ and $Q_{s}\tilde{\tau}_{0}=1$. For all initial distributions, the longitudinal hard scale shows an identical classical scaling behavior.

III.2 Effects of the time dependence of the Coulomb logarithm

In some of the earlier studies that employ the kinetic equation in the small angle approximation, the Coulomb logarithm (14) is taken to be a time independent constant Blaizot:2013lga ; Blaizot:2014jna . In the present study, since the Debye mass given by Eq. (15) is a time dependent IR scale, the Coulomb logarithm varies in time. We will investigate the effect of this time dependence in this subsection.

This is demonstrated in Fig. 8 for the parameters $g=10^{-2}$, $n_{0}=0.1$, $\xi_{0}=2$ and $Q_{s}\tilde{\tau}_{0}=1$. As the Debye mass decreases in time as $m_{D}^{2}\sim 1/\tau$, the Coulomb logarithm increases in time logarithmically. The increase is slow except for early times. Therefore the classical scaling behavior is not spoiled by the time dependence of the IR cutoff. The only effect is the change of the effective interaction strength. Because $\mathcal{L}$ appears as an overall multiplicative factor in the collision terms, the increase of $\mathcal{L}$ indicates that the time evolution is faster. This can be confirmed in Fig. 9, where the longitudinal hard scale is plotted for the case where the time dependence of $\mathcal{L}$ is taken into account and compared to results where $\mathcal{L}$ is fixed to its initial value555In the computation of the latter case, we have used $N_{\kappa}=512$ for the number of angular grid points in order to precisely calculate $\Lambda_{L}^{2}$ at later times.. The transition from the free streaming regime to the classical scaling regime happens later in the constant $\mathcal{L}$ computation, indicating the effective interaction strength is weaker in this case.

IV Numerical results II. $N_{f}=3$

Now that we have confirmed that the kinetic equation within the small angle approximation can qualitatively describe the classical scaling behavior of overpopulated gauge fields, we will add $N_{f}=3$ quark degrees of freedom to the expanding Glasma. We will assume that the initial distribution for quarks is given by

[TABLE]

The shape of this initial distribution is irrelevant for the following discussion on the classical scaling behavior; for simplicity, we have therefore chosen the same momentum dependence as that for the initial gluon distribution. We note that on account of the Pauli principle, the initial occupancy $F_{0}$ cannot exceed unity.

IV.1 Self-similar evolution

Figure 10 displays the time evolution of the quark number density multiplied by time. Three different values for the initial quark occupancies $F_{0}$ are compared. The parameters used in these computations are $g=10^{-2}$, $n_{0}=0.1$, $\xi_{0}=2$, and $Q_{s}\tilde{\tau}_{0}=1$. If quark number is conserved, $\tilde{\tau}n_{q}$ should stay constant. The increase in this quantity with time, seen in Fig. 10, clearly indicates the effect due to quark pair production via gluon fusion. The contribution from this process is however rather mild because the effective interaction strength with gluons is governed by the factor $g^{2}f_{0}=n_{0}$; as our choice of $n_{0}$ indicates, this quantity is small in the regime where the kinetic description is valid. On the other hand, at early times in heavy-ion collisions, \tau\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}Q_{s}^{-1}, the factor $g^{2}f_{0}$ is of order unity and quark pair production can be quite large due to nonperturbative processes, as hinted at by extant real-time lattice simulations Gelis:2005pb ; Gelfand:2016prm ; Tanji:2016dka .

The second moment of the quark distribution evaluated at $p_{z}=0$ is plotted in the left panel of Fig. 11 for different times. The initial quark occupancy is chosen to be $F_{0}=0.5$, and the other parameters are the same as those for Fig. 10. Here we show only the quark distribution because the gluon distribution is little affected by the existence of quarks for Glasma initial conditions. We observe that, just as in the gluon case, the typical transverse momentum is centered around $Q_{s}$. In the right panel of Fig. 11, we plot the quark distribution multiplied by the factor $(Q_{s}\tilde{\tau})^{2/3}$. At later times, Q_{s}\tilde{\tau}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}20, the rescaled distributions at different times overlap nicely. This indicates that the quark transverse momentum distribution obeys the same scaling law as the gluon transverse momentum distribution.

It is not just the transverse momentum distributions that satisfy scaling laws. The quark longitudinal momentum distributions satisfy the same scaling law as that for gluons. This is displayed in Fig. 12. These observations can be summarized into the scaling expression,

[TABLE]

which has the identical form as that of Eq. (50) for gluons. This scaling behavior is further confirmed in Fig. 13, where the time evolution of the longitudinal hard scale is plotted for quark distributions with different values of the initial occupancy $F_{0}$. We observe that a transition from free streaming behavior to the scaling form occurs for any value of $F_{0}$.

At first sight, the exhibition of scaling behavior by the quark distribution might appear strange because Bose enhancement is a prerequisite for the emergence of the scaling behavior. However this is in fact not at all strange because what needs to be Bose-enhanced is a scattering amplitude. When quarks elastically scatter with gluons, as is transparent from the kinetic equations, the presence of gluons in the final state causes the process to be Bose-enhanced. In both the equation for gluons (19), and that for quarks (20), the diffusion terms involve the common factor $I_{a}$, which can be interpreted as the density of scatterers. The integral $I_{a}$ in the Glasma is dominated by the gluon distribution, and it is therefore enhanced by the Bose factor. Thus quarks undergo small angle scattering processes that are Bose-enhanced just as for gluons, and therefore exhibit the same scaling behavior.

To demonstrate the emergence of scaling behavior for the quark distributions more explicitly, we will follow the argument for gluons presented in Ref. Berges:2013fga . In the classical scaling regime of the Glasma, the typical longitudinal momentum is much smaller than the typical transverse momentum. Therefore, the collision term in the small angle approximation is dominated by the derivative term in $p_{z}$. Furthermore, as observed in Fig. 10, quark pair production is slow in the kinetic regime, and the source term may be neglected. The kinetic equation for quarks (Eq. (20)) can therefore be approximated by

[TABLE]

where we define $\hat{q}$ to be

[TABLE]

a quantity proportional to the integral $I_{a}$ and dominated by the gluon distribution. Since the gluon distribution has the scaling form specified in Eq. (50), the time dependence for $\hat{q}$ can be identified as

[TABLE]

up to the logarithmic time dependence contained in $\mathcal{L}$. If we assume that the following scaling form of the quark distribution

[TABLE]

is a stationary solution of the simplified kinetic equation (56), it follows that

[TABLE]

Furthermore, from the approximate number conservation for quarks, $n_{q}\sim 1/\tau$, one obtains a constraint

[TABLE]

Since the typical longitudinal momentum is much smaller than the typical transverse momentum, the longitudinal pressure density $\mathcal{P}_{L}$ is negligible compared with the energy density $\mathcal{E}$. In this case, the energy density approximately behaves as $\mathcal{E}\sim 1/\tau$, from which it follows that

[TABLE]

By these, the remaining exponents are determined to be

[TABLE]

To summarize, quark distributions show the same scaling behavior as that exhibited by the gluon distribution when (i) the diffusion constant $I_{a}$ is dominated by overoccupied gluons, (ii) the typical longitudinal momentum is much smaller than the typical transverse momentum, and (iii) the quark number density $\tau n_{q}$ does not vary rapidly in time.

Though quark distributions qualitatively obey the same scaling behavior as those for gluons, there are quantitative differences in their temporal evolution. In Fig. 14, the longitudinal hard scale is compared for gluons and quarks ($F_{0}=0.5$). The transition from the free streaming behavior to the classical scaling behavior happens later for quarks than gluons. This can be understood simply as a consequence of the smaller effective coupling of quarks seen in the comparison between Eqs. (19) and (20). A similar tendency is observed in the bulk anisotropy (the ratio of the longitudinal pressure to the transverse pressure) plotted in Fig. 15. The quark distribution is more anisotropic relative to the gluon distribution as quarks scatter more weakly than gluons.

IV.2 Approach to equilibrium

In the bottom-up thermalization scenario, the inelastic collinear splitting process that produces soft gluons plays a crucial role in thermalization. As the kinetic equations (19) and (20) contain only elastic collision terms, they are applicable only to the classical scaling regime in bottom-up thermalization.

Nevertheless, albeit not justified on parametric grounds, it is an amusing exercise to crank up the coupling and see how both kinetic and chemical equilibration occur via only elastic collisions. This is for instance the strategy adopted in a number of phenomenological models of heavy-ion collisions – and it is useful to understand, from this perspective at least, how equilibration occurs. With only elastic collisions included , thermalization takes a parametrically long time $\sim Q_{s}^{-1}\exp(\alpha_{s}^{-1/2})$ Mueller:1999fp ; Bjoraker:2000cf ; Serreau:2001xq – much longer than that in the bottom-up scenario.

In order to observe the approach to equilibrium in a practical computational time, we take a relatively large value for the coupling $g=0.5$ and the gluon occupation coefficient $n_{0}=3/4$, where admittedly the validity of the kinetic equations is marginal. Furthermore, we use the initial anisotropy parameter $\xi_{0}=1$ and a relatively large initial time $Q_{s}\tilde{\tau}_{0}=100$.666These conditions correspond to the initial time $Q_{s}\tau_{0}\simeq 100$ in terms of the original time variable. For the quark initial occupation number, we compare different values $F_{0}=0$, 0.5 and 1. When the system gets close to equilibrium, the assumption of a constant UV cutoff $q_{\text{max}}=Q_{s}$ in the Coulomb logarithm is inadequate; we will instead use the mean transverse momentum in Eq. (16) as the cutoff.

The ratio of pressures, $\mathcal{P}_{L}/\mathcal{P}_{T}$, is plotted as a function of time in Fig. 16. The longitudinal pressure $\mathcal{P}_{L}$ and the transverse pressure $\mathcal{P}_{T}$ include here the contributions from both gluons and quarks. The ratio decreases at early times due to the longitudinal expansion–unlike the result for weaker couplings however, already at $g=0.5$ one observes a turnaround in this ratio and the ratio begins to increase. While full isotropization is not achieved777It should be noted though that the value of the coupling, $g=0.5$, that we choose, $Q_{s}$ is enormous and is far above the TeV scale – hence values of $Q_{s}\tau\sim 1000$ on plots still correspond to times far less than a Fermi. The fact that the magnitude $Q_{s}$ should be commensurate with the value of $g$ chosen is often obscured in a number of works, where $Q_{s}\sim 1$ GeV is assumed even for very weak couplings. by $Q_{s}(\tilde{\tau}-\tilde{\tau}_{0})=3000$, the ratio reaches 0.8.

In Fig. 17, the ratio of the quark number and the gluon number is plotted as a function of time. In thermal equilibrium with vanishing chemical potential, the ratio should be $3N_{c}N_{f}/(2(N_{c}^{2}-1))=27/16$, which is shown as a gray dashed line. The ratio of quark and gluon densities approaches its equilibrium value more quickly relative to the ratio of pressures. The simple message is that while pressure isotropization requires the system to overcome the effect of expansion, in contrast, equipartition is unhindered since both quarks and gluons experience the expansion.

V Summary and outlook

We presented in this paper numerical results from the solution of kinetic equations for an expanding Glasma of overoccupied gluons and quarks. We considered only elastic 2$\leftrightarrow$2 scatterings amongst the quarks and gluons. As may be anticipated, we observed that a cascade of particle number to the infrared develops generating an enhancement in this region; this phenomenon can be interpreted heuristically as signaling the formation of a transient Bose-Einstein condensate. For $N_{f}=0$, and for small values of the coupling where classical assumptions are robust, we find that our simple kinetic equation for the Glasma qualitatively reproduces key features of classical-statistical real-time lattice simulations of the 3+1-dimensional Yang-Mills equations. These numerical simulations, for sufficiently small couplings, should completely capture the physics of the expanding Glasma.

The fact that results of our kinetic simulations reproduce the key features of the classical-statistical simulation results is truly remarkable. This is because, on the face of it, these kinetic equations do not include $2\leftrightarrow 3$ scatterings that are also present in the classical stage of bottom-up thermalization. Such inelastic scatterings would, in principle, prevent a BEC from developing Kurkela:2011ti ; York:2014wja ; Blaizot:2016iir . Since a first principles weak coupling kinetic treatment should contain these contributions, the bottom-up kinetic scenario does not include the presence of a transient BEC. However the bottom-up framework is also not fully robust because, as we noted previously, it does not contain the effects of late time plasma instabilities that would modify the temporal evolution of the Glasma Kurkela:2011ti . As discussed in detail in Berges:2015ixa , the fact that the bottom-up scenario captured some of the key features of the numerical simulations posed a challenge for conventional ab initio kinetic treatments of the Glasma.

Further, as also noted in Berges:2015ixa , the bottom-up framework was unable to reproduce the temporal behavior of the anisotropy measure $P_{L}/P_{T}$ seen in the numerical simulations. From the bottom-up based simulations presented in the literature Kurkela:2015qoa , it is unclear whether the $1/p_{\!\perp}$ dependence of the single particle distributions is reproduced. The fact that our simplified kinetic treatment with only 2$\leftrightarrow$2 scattering captures qualitative features of the overoccupied Glasma therefore confirms the conjecture in Berges:2015ixa that the infrared dynamics is playing a bigger role than anticipated not only in suppressing late time plasma instabilities but in the evolution of the bulk anisotropy as well.

The possibility that a transient BEC might provide an effective description of infrared dynamics in the Glasma was conjectured in Blaizot:2011xf and developed further in Blaizot:2013lga . However these estimates were made without the benefit of detailed kinetic simulations of the expanding Glasma and assumed a constant value of $P_{L}/P_{T}$ even for weak couplings. This constant value for the anisotropy was disfavored by classical-statistical simulations. It is therefore perhaps ironic that an effective BEC description of the classical stage of the Glasma survives in a qualitative comparison of our effective kinetic description and the full classical-statistical Yang-Mills simulations. It will be important in future to quantify this comparison and understand better how BEC formation is influenced by $2\leftrightarrow 3$ processes at weak coupling.

We note that for self-interacting scalar theories, the existence of a BEC is firmly established, for the same expanding geometry as the Glasma Berges:2014bba . A big puzzle discussed in Berges:2014bba was why the simulations of scalar theories showed identical universal behavior to that of the Glasma. A conventional kinetic description of a scalar theory would not capture its dynamics. In contrast, an effective kinetic description with 2$\leftrightarrow$2 scattering, albeit with couplings modified by infrared dynamics, reproduces the infrared behavior in numerical simulations in a controlled fashion Orioli:2015dxa . This is similar to our BEC+2$\leftrightarrow$2 scattering effective description albeit it should be noted that our BEC will not capture important features of the IR dynamics in the scalar case. It is challenging to perform controlled computations of infrared dynamics in a gauge theory as have been performed for a scalar theory, or even interpret this IR dynamics definitively as a BEC. However, as we noted previously, recent computations of the spatial string tension in the Glasma Mace:2016svc ; Berges:2017igc provide an essential guidepost for future studies.

The other key result of this paper was the first numerical simulation of the temporal evolution of quark distributions in the expanding overoccupied Glasma. We showed that light quark flavors obey identical scaling distributions to those of gluons. Quarks however take longer to achieve the scaling behavior because their effective interactions are weaker than those of gluons. The scaling results have already been employed to estimate the photon yield in the Glasma Berges:2017eom ; our results there showed that photon yields from the Glasma are significant relative to those from the thermalized QGP. We also studied the rate at which chemical equilibration is achieved. This result is sensitive to the initial quark occupancies. For maximal initial quark occupancies, chemical equilibration is approached rapidly; it remains to be seen how inelastic collisions which become influential for times \tau\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}Q_{s}^{-1}\,\alpha_{s}^{-3/2} modify this result.

As an outlook for future work, it will be important to perform more detailed comparisons of our effective kinetic framework both with classical-statistical Yang-Mills simulations and kinetic simulations in the bottom-up framework that include inelastic scattering contributions. These comparisons are relevant both for our results for $N_{f}=0$ and at $N_{f}=3$. Further, a deeper theoretical understanding of our results, in the context of the puzzles outlined is desirable. A useful guide is the body of work employing two-particle irreducible $1/N$ techniques that have proved powerful in the context of scalar theories Berges:2015kfa . These also provide guidance in the extrapolation of weak coupling techniques to couplings relevant for heavy-ion collisions Berges:2016nru . Finally, we note that in our discussion of photon yields in Berges:2017eom , we did not consider the contribution to photon production from the interaction of quarks with the BEC. Such a possibility was discussed previously in Chiu:2012ij . It will be useful to reconsider the possible effects of a transient BEC in this context, and in other contexts Begun:2015ifa ; Gangadharan:2016vez ; Osada:2017oxe , such as for instance the intriguing results on Hanbury-Brown–Twiss correlations from the ALICE experiment at the LHC Adam:2015pbc .

Acknowledgements

We would like to thank Juergen Berges, Kirill Boguslavski, Aleksi Kurkela and Alexander Rothkopf for valuable discussions and comments. R. V. is supported under DOE Contract No. DE-SC0012704. He would like to thank the Institut für Theoretische Physik, Universität Heidelberg for kind hospitality and support via the Excellence Initiative during the early stages of this work. This work is part of and supported by the DFG Collaborative Research Centre “SFB 1225 (ISOQUANT)”.

Appendix A Numerical method

In this appendix, we will discuss details of the numerical method employed to solve the kinetic equations. Since the kinetic equation for gluons Eq. (19), and that for quarks Eq. (20), have the same structure, it is sufficient to discuss only Eq. (19). In terms of the variables $p=\sqrt{p_{\!\perp}^{2}+p_{z}^{2}}$ and $\kappa=p_{z}/p$, Eq. (19) can be rewritten as

[TABLE]

where $\mathcal{R}$ denotes the (rescaled) flow in the radial direction,

[TABLE]

The form of this flow is exactly the same as that in earlier studies for isotropic systems, except that now $\mathcal{R}$ also depends on the angular variable $\kappa$. The second term on the r.h.s of Eq. (A) describes diffusion in $\kappa$–this term is of course absent in the isotropic case.

To solve the kinetic equations numerically, we discretize the momentum variables $p$ and $\kappa$ into $N_{p}$ and $N_{\kappa}$ grids, respectively, as

[TABLE]

and

[TABLE]

Before discretizing Eq. (A), we rewrite $p\partial\tilde{f}/\partial p$ that appears on the l.h.s as

[TABLE]

This “integration by parts” is helpful in preserving number conservation (Eq. (33)) to high numerical accuracy. We discretize Eq. (A) as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

In Eqs. (72) and (76), $p_{0}$ should be set to zero. For this choice in the discretization for the equation of motion, the integrals $\tilde{I}_{a}$ and $\tilde{I}_{b}$ can be replaced by

[TABLE]

where $\sum_{j,\text{trap}}$ denotes the trapezoidal sum,

[TABLE]

Note that the discretization of $\int\!dp$ is different for $\tilde{I}_{a}$ and $\tilde{I}_{b}$. This is done in order to preserve energy conservation (Eq. (34)) to high accuracy. For the discretization of $\tilde{I}_{c}$, there is no guidance from conservation laws. For simplicity, we employ the same discretization as for $\tilde{I}_{b}$:

[TABLE]

As noted, the kinetic equation for quarks, Eq. (20), is discretized exactly in the same way.

To solve Eq. (A), a boundary condition for the radial flow $\mathcal{R}_{ij}$ at $i=0$ is necessary. Since we consider only number conserving elastic collisions, the onset of the Bose-Einstein condensate occurs at a finite time for an overoccupied initial condition Blaizot:2013lga . When the condensate is generated, special care is necessary for the boundary condition at $p=0$. We adopt the method presented in Ref. Blaizot:2014jna . Before the onset of the BEC, the gluon distribution is less singular than $1/p$. Therefore, $\lim_{p\to 0}\mathcal{R}(p,\kappa)=0$. Once the BEC is formed, the distribution shows a $1/p$ behavior in the infrared, and the flow into the deep infrared becomes nonzero:

[TABLE]

where $c$ is the coefficient of the $1/p$ term of the gluon distribution; $\tilde{f}(p)\simeq c/p$. In numerical computations, we employ the following boundary condition

[TABLE]

where the coefficient $c_{j}$ is determined by

[TABLE]

and $\tau_{c}$ is the time for the onset of condensation, defined as the instant that $\tilde{I}_{b}c_{j}^{2}$ becomes larger than $\tilde{I}_{a}c_{j}$. In principle, $\tau_{c}$ can depend on $j$. However for sufficiently small $p_{\text{min}}$, it does not depend on $j$. For quarks, the radial flow at the origin is always vanishing because the distribution for fermions is bounded to be below unity. Therefore the boundary condition for the quark kinetic equation is $\mathcal{R}_{0j}=0$.

For the time evolution, we employ the alternative direction implicit method. However, we treat the nonlinear terms in $f$ and $F$ explicitly.

Appendix B Bose-Einstein condensate

The Landau kinetic equation for the 2$\leftrightarrow$2 elastic collisions was employed previously to discuss the onset of Bose-Einstein condensation Blaizot:2013lga ; Blaizot:2014jna . However number changing inelastic processes play a crucial role in the infrared region and they may hinder the onset of the condensate Kurkela:2011ti ; York:2014wja ; Blaizot:2016iir . While they will ensure that any BEC formed is transient in the “long run”, a BEC or more generally a large IR occupancy may play a role in the classical regime that is the focus of this study. For the expanding Glasma, a role of inelastic processes in BEC formation needs more careful study, for a range of couplings, but is outside the scope of the investigation here.

Here we will focus on the numerical implementation of BEC formation in 2$\leftrightarrow$2 scattering and discuss its properties. Although an accurate description of condensate formation and IR dynamics is outside of the range of applicability of our kinetic equations, BEC formation affects the applicability of our kinetic equations because of the boundary condition in momentum space, as discussed in the previous section. For simplicity, we will show only results for $N_{f}=0$ – without quarks. As long as gluons are initially overoccupied, the effect of quarks on the gluon condensate is negligible. The effects of quarks for moderate values of $f_{0}$ in a non-expanding system have been discussed in Ref. Blaizot:2014jna .

As discussed in the previous section, the condensate forms when the flow to the deep infrared (85) sets in. In Fig. 18, we plot the flow, averaged over angle in momentum space,

[TABLE]

as a function of time. $c_{j}$ is computed by (86). Results with different $p_{\text{min}}$ are compared888For $p_{\text{min}}/Q_{s}=10^{-2}$, we have used $N_{p}=500$. For other values of $p_{\text{min}}$, we adjusted $N_{p}$ so that $\Delta u$ is the same.. At early times, we observe a positive flow corresponding to a UV cascade from infrared to higher momenta. This is an artificial result that occurs because we have nonzero $p_{\text{min}}$. Indeed, for smaller $p_{\text{min}}$, positive flow is suppressed. This observation justifies the prescription (85).

After a time $Q_{s}\tilde{\tau}_{c}\simeq 2$, the flow becomes negative, indicating particle flow into the infrared region $p<p_{\text{min}}$. For sufficiently small $p_{\text{min}}$ (p_{\text{min}}/Q_{s}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}10^{-2}), the time $\tilde{\tau}_{c}$ and the behavior of $\bar{\mathcal{R}}_{0}$ after $\tilde{\tau}=\tilde{\tau}_{c}$ are insensitive to the value of $p_{\text{min}}$.

Figure 19 shows the time evolution of the BEC number density, defined as

[TABLE]

Curves in the plot are insensitive to the value of $p_{\text{min}}$ for sufficiently small $p_{\text{min}}$. Since $\tilde{\tau}\tilde{n}_{\text{cond}}$ can be regarded as the total particle number in the momentum sphere $p<p_{\text{min}}$, the insensitively to $p_{\text{min}}$ indicates that these particles condense to a zero mode.

The time evolution of particle number in the BEC is plotted in Fig. 20 for different values of the coupling. The initial occupancy coefficient is fixed to be $n_{0}=0.1$. As we parametrize the initial occupancy to be $f_{0}=n_{0}/g^{2}$, different values of the coupling correspond to different initial occupation number. Other parameters used in these calculations are $\xi_{0}=1$ and $Q_{s}\tilde{\tau}_{0}=1$. In the figure, $\tilde{\tau}n$ divided by $f_{0}$ is plotted. For $g=0.5$, which corresponds to $f_{0}=0.4$, $\tau n_{\text{cond}}$ stays zero indicating that the condensation does not occur. For weaker coupling (larger initial occupancy), condensation occurs around $Q_{s}(\tilde{\tau}-\tilde{\tau}_{0})=1$ and $\tau n_{\text{cond}}$ begins increasing at earlier times. In a nonexpanding system, the number of condensate particles approaches its equilibrium value, which can be nonzero for overoccupied initial conditions999One expects of course that quantum corrections at late times will destroy the BEC.. In contrast, condensation disappears in the expanding system, simply by its dilution from the expansion. Indeed, the increase in the condensate density stops at some point and a turnover is seen for $g=0.05$ and 0.1.

We plot in Fig. 21 the distribution function before and after the onset of the condensate. The transverse momentum distribution evaluated at $p_{z}=0$ and the longitudinal momentum distribution at $p_{\!\perp}=0$ are compared. The initial conditions $g=10^{-2}$, $n_{0}=0.1$ ($f_{0}=10^{3}$), $\xi_{0}=1$ and $Q_{s}\tilde{\tau}_{0}=1$ correspond to an overoccupied and isotropic distribution. The time for the onset of condensation is $Q_{s}\tilde{\tau}_{c}\simeq 2.08$. Although the distribution is initially isotropic, the hard sector p\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle>}}{{\sim}}>}Q_{s} becomes quickly anisotropic due to the longitudinal expansion of the system. On the other hand, the soft sector p\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}0.1Q_{s} stays isotropic once the condensate is formed and maintains a $p^{-1}$ shape after $\tilde{\tau}=\tilde{\tau}_{c}$. However the region where the $p^{-1}$ shape is realized shrinks in the longitudinal direction due to the expansion.