TThe First fm/c of
Heavy-Ion Collisions
S. Schlichting, D. Teaney, Fakult¨at f¨ur Physik, Universit¨at Bielefeld, D-33615 Bielefeld, Germany; email:[email protected] Department of Physics and Astronomy, Stony Brook University, Stony Brook,New York 11794, USA; email: [email protected]. Xxx. Xxx. Xxx. 2019. AA:1–31https://doi.org/10.1146/((please addarticle doi))Copyright c (cid:13)
Keywords
Heavy-Ion Collisions; Quark-Gluon Plasma; QCD Kinetic Theory;Thermalization;
Abstract
We present an introductory review of the early time dynamics of high-energy heavy-ion collisions and the kinetics of high temperature QCD.The equilibration mechanisms in the quark-gluon plasma uniquely re-flect the non-abelian and ultra-relativistic character of the many bodysystem. Starting with a brief expose of the key theoretical and ex-perimental questions, we provide an overview of the theoretical toolsemployed in weak coupling studies of the early time non-equilibriumdynamics. We highlight theoretical progress in understanding differ-ent thermalization mechanisms in weakly coupled non-abelian plasmas,and discuss their relevance in describing the approach to local thermalequilibrium during the first fm /c of a heavy-ion collision. Some impor-tant connections to the phenomenology of heavy-ion collisions are alsobriefly discussed. a r X i v : . [ nu c l - t h ] A ug ontents
1. Introduction and motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22. Early time dynamics of heavy-ion collisions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.1. Microscopics of the initial state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2. Bottom-up equilibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63. QCD Kinetics: a brief review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73.1. Elastic scattering and momentum diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83.2. Collinear radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104. Basics of weak coupling thermalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134.1. Overoccupied systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134.2. Underoccupied systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164.3. Generalization to anisotropic systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205. Simulations of early time dynamics and heavy-ion phenomenology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235.1. Approach to hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235.2. Quark production and chemical equilibration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255.3. Small scale fluctuations and pre-flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 266. Outlook and small systems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
1. Introduction and motivation
The purpose of ultra-relativistic heavy ion collisions is to produce and to characterize theproperties of the Quark-Gluon-Plasma (QGP), which is an extreme state of Quantum-Chromo-Dynamic (QCD) matter that was also present in the early universe, during thefirst mircoseconds after the big bang. Over the last two decades, experiments at the Rel-ativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) have collideda variety of nuclei over a wide range of energies, and, at least in the collisions of largenuclei, these experiments show that the produced constituents re-interact, and exhibitmulti-particle correlations with wavelengths which are long compared to the microscopiccorrelation lengths, providing overwhelming evidence of collective hydrodynamic flow (1).Hydrodynamic simulations of these large nuclear systems describe the observed correlationsin exquisite detail with a minimal number of parameters (1). In smaller systems such asproton-proton (pp) and proton-nucleus (pA) long range flow-like correlations amongst theproduced particles have also been observed (2, 3), and these observations drive current re-search into the equilibration mechanism of the QGP. This research aims to understand howthe observed correlations change with system size, and approach the hydrodynamic regimefor large nuclei.Explaining approximately how an equilibrated state of quarks and gluons emerges fromthe initial wave functions of the incoming nuclei has been one of the central goals of the heavyion theory community for a long time. Even though genuinely non-perturbative real-timeQCD calculations are currently not available to address this question (as they suffer from asevere sign problem), significant progress has been achieved in understanding properties ofthe initial state and the equilibration mechanism based on ab-initio calculations at weak andstrong coupling. Here we focus on the weak coupling description, based on the idea that athigh energy density and temperatures the coupling constant between quarks and gluons α s becomes small, and weak coupling methods can be used to analyze the initial production ofquarks and gluons, and the kinetic processes which ultimately lead to a thermalized QGP. hen extrapolated to realistic coupling strength, the weak coupling approach based onperturbative QCD and strong-coupling approaches based on the holography yield similarresults for the macroscopic evolution of the system (5). For a recent review of the strongcoupling description we refer to (4).The weakly coupled picture of the equilibration process in high energy collisions wasoutlined in a seminal paper by Baier, Mueller, Schiff and Son (BMSS) (6), and is referredto as the bottom-up thermalization scenario,which is schematically depicted in Fig. 2. Weprovide a short review of bottom-up in Sect. 2, and then describe recent reanalyses whichhave clarified and extended the original picture considerably. These extensions have turnedthe parametric estimates of BMSS into hard numbers, which can be used to make contactwith the experimental data.We emphasize that the study of the equilibration mechanisms in non-abelian gaugetheories, such as QCD, is of profound theoretical interest, and much of the research intothermalization is only tangentially driven by the immediate needs of experimental heavyion physics program. In this spirit this review aims to cover some of the most importanttheoretical developments regarding the equilibration mechanism in non-abelian plasmas.Starting with an introductory discussion of the basic physics picture of the early stages ofhigh-energy heavy-ion collisions in Sect. 2, the subsequent sections, Sects. 3 and 4, provide amore detailed theoretical discussion of the underlying theory and the equilibration process ofweakly coupled non-abelian plasmas. New developments based on microscopic simulationsand connections to heavy-ion phenomenology are then discussed in Sect. 5.
2. Early time dynamics of heavy-ion collisions
When two nuclei collide at high energies, they pass through each other scarcely stopped,leaving behind a debris of highly excited matter which continues to expand longitudi-nally (7). Since the system is approximately invariant under boosts in the longitudinal( z ) direction, one point functions of the stress tensor and other fields in the central rapidityregion, i.e. the region close to the original interaction point, only depend on proper time τ = √ t − z , but do not depend on the space time rapidity η = log(( t + z ) / ( t − z )). Inco-moving ( τ, x, y, η ) coordinates, the metric is ds = − dτ + dx + dy + τ dη , (1)indicating that the boost-invariant system is continually expanding along the beam axis, dz = τ dη . While initially the system is in a state far from local thermal equilibrium,phenomenology suggests that on a time τ hydro ∼ / c the plasma of quarks and gluons issufficiently close to equilibrium that hydrodynamic constitutive relations are approximatelysatisfied and the subsequent evolution can be described with hydrodynamics.While the longitudinal structure is approximately homogenous in space time rapidity,the transverse structure of the fireball is always inhomogenous, reflecting the initial geom-etry of the collision. In Fig. 1 we show a typical transverse (entropy density) profile that isused to initialize hydrodynamic simulations of the space time evolution. While the averagegeometry is characterized by the nuclear radius R A , one finds that in any realistic event-by-event simulation there are smaller length scales in the initial geometry of order the protonradius, R p (cid:28) R A , which arise from fluctuations in the positions of the incoming protons.Such geometric fluctuations are responsible for many of the most prominent flow observ-ables in heavy ion collisions such as e.g. the triangular flow (8). Still smaller fluctuations ••
When two nuclei collide at high energies, they pass through each other scarcely stopped,leaving behind a debris of highly excited matter which continues to expand longitudi-nally (7). Since the system is approximately invariant under boosts in the longitudinal( z ) direction, one point functions of the stress tensor and other fields in the central rapidityregion, i.e. the region close to the original interaction point, only depend on proper time τ = √ t − z , but do not depend on the space time rapidity η = log(( t + z ) / ( t − z )). Inco-moving ( τ, x, y, η ) coordinates, the metric is ds = − dτ + dx + dy + τ dη , (1)indicating that the boost-invariant system is continually expanding along the beam axis, dz = τ dη . While initially the system is in a state far from local thermal equilibrium,phenomenology suggests that on a time τ hydro ∼ / c the plasma of quarks and gluons issufficiently close to equilibrium that hydrodynamic constitutive relations are approximatelysatisfied and the subsequent evolution can be described with hydrodynamics.While the longitudinal structure is approximately homogenous in space time rapidity,the transverse structure of the fireball is always inhomogenous, reflecting the initial geom-etry of the collision. In Fig. 1 we show a typical transverse (entropy density) profile that isused to initialize hydrodynamic simulations of the space time evolution. While the averagegeometry is characterized by the nuclear radius R A , one finds that in any realistic event-by-event simulation there are smaller length scales in the initial geometry of order the protonradius, R p (cid:28) R A , which arise from fluctuations in the positions of the incoming protons.Such geometric fluctuations are responsible for many of the most prominent flow observ-ables in heavy ion collisions such as e.g. the triangular flow (8). Still smaller fluctuations •• The First fm/c of Heavy-Ion Collisions 3 x R p R A /Q s c ⌧ hydro p τ hydro s τ h y d r o s ( τ h y d r o ) ( G e V ) x (fm)one eventaverage Figure 1 (left) Illustration of the two nuclei as they are passing through each other. Classical color fieldconfigurations just after passage were described in (9, 10, 11) and feature strong longitudinalchromo-electric and chromo-magnetic fields, which rapidly decohere on a timescale of ∼ Q s (12).(right) Snapshot of a typical entropy density profile used in event-by-event hydrodynamicsimulations of heavy-ion collisions (13). Smaller scale fluctuations on microscopic length scales ∼ /Q s are not shown, and are indicated by the black dot. of order the inverse saturation momentum Q − s (see Sect. 2.1) are not shown in this figure.Different scales in Fig. 1 should be compared to the distance scale cτ hydro , which providesan estimate of the causal propagation distance during the approach to equilibrium. Wewill generally assume that cτ hydro is short compared to the nuclear radius, cτ hydro (cid:28) R A ,such that on average the transition from the non-equilibrium state towards thermal equi-librium proceeds locally in space and can discussed at the level of individual cells of size cτ hydro . Short distance fluctuations on scales cτ hydro spoil this picture; however such effectswere neglected in the original bottom-up scenario and we will follow this assumption byapproximating the evolution of the system as homogenous in transverse space and spacetime rapidity throughout most of this review. Shortcomings of this approximation will bediscussed further in Sec. 5 and 6 along with recent extensions of the original work of BMSS,which incorporate short distance fluctuations of the nucleon positions on scales R p ∼ cτ hydro into the description of the first fm/c of heavy-ion collisions. In each small circle of size cτ hydro in the transverse plane the initial production of quarks andgluons in momentum space follows from the Color-Glass-Condensate (CGC) effective theoryof parton saturation (14, 15). Briefly, in this theory the incoming nuclei are highly lengthcontracted by an ultra-relativistic factor γ (cid:29)
1, and the density of gluons per transversearea and rapidity in the wave functions of the nuclei, ( dN/dy ) /πR A , grows with increasingcollision energy. Here dN/dy is the number of gluons per rapidity y which is related toBjorken x bj , dy = dx bj /x bj . This transverse density of gluons determines a momentumscale, known as the saturation momentum Q s , which at very high energies can becomelarge compared to Λ QCD Q s ∼ α s πR A dNdy (cid:29) Λ QCD . (2) he saturation momentum Q s sets the momentum scale for the transverse momentumdistribution of partons in the wave functions. For Q s (cid:29) Λ QCD the coupling constant issmall α s ( Q s ) (cid:28)
1, and the evolution of the system can be treated using weakly coupledmethods. Further, the number of gluons per phase space cell in the incoming wave functionsis large 1 πQ s R A dNdy ∼ α s (cid:29) , (3)and in this regime the evolution of the system is classical. Thus, the production of gluonsand their initial evolution system is determined by solving the non-linear classical Yang-Mills equations of motion (9, 10, 11). In practice, the saturation momentum is Q s ∼ QCD there willalways be important quantum corrections to the CGC formalism, which will almost becompletely neglected in this review.In an important set of papers, the initial conditions for the classical fields in the forwardlight cone just after the intitial crossing of the two nuclei were worked out (analytically) bymatching the classical fields just before the collision with those just after crossing (9, 10, 11).These initial conditions consist of strong longitudinal fields, E z and B z , which as illustratedin Fig. 1 is somewhat reminiscent of a parallel plate capacitor (12). Indeed, the averagestress tensor for a boost invariant, or Bjorken, expansion and a conformal system (with T µµ = 0) must take the form (cid:104) T µν (cid:105) = ( − (cid:15), P T , P T , P L ) , with (cid:15) = 2 P T + P L . The matchingprocedure (9, 10, 11) shows that (cid:104) T µν (cid:105) = ( − (cid:15), (cid:15), (cid:15), − (cid:15) ), and thus, the initial longitudinal“pressure” P L is negative as for a constant electric (or magnetic) field in the z -direction inclassical electrodynamics. These strong longitudinal fields rapidly decrease on a time scaleof ∼ Q s as the classical field configuration decoheres.The initial conditions outlined in the preceding paragraph motivated the first classicalsimulations of gluon production in the longitudinally expanding boost invariant geome-try (16, 17). In the original formulations the classical fields were assumed to remain effec-tively 2+1 dimensional, i.e. strictly independent of rapidity as a function of time τ , reflectingthe fact that the initial conditions are boost invariant up to quantum corrections of order α s . However, such quantum fluctuations provide the seed from which classical instabilitiesdevelop in the longitudinal direction (18, 8), such that the gluonic fields quickly becomechaotic in all three dimensions and the classical solutions are only rapidity-independent onaverage. The instabilities grow as e Γ √ Q s τ , with Γ ∼
1, limiting the applicability of strictlyboost invariant simulations to short times, τ (cid:46) Q − s log (1 /α s ) (19, 20, 21). In spite ofthis shortcoming, strictly boost invariant simulations of classical field dynamics form thebasis of phenomenological studies of particle production and early time dynamics in theIP-Glasma model (22, 23).During the classical evolution the field strength decreases due to the longitudinal ex-pansion, and eventually the equations of motion linearize. For times long enough τ Q s (cid:29) α − s ( Q s ). In this regime, either kinetic theory orclassical field theory can be used to simulate the evolution of the system (24, 25, 26). Inparticular, it is sensible to talk about the gluon phase space distribution, as opposed to theclassical field configuration. The initial phase space distribution of gluons f ( τ, x , p ) canbe determined from the classical simulations by evaluating the Wigner transform of equaltime two point functions of gauge fields, after fixing a physical gauge such as the CoulombGauge (see for instance ref. (27, 28, 29)). Due to the longitudinal expansion of the system, ••
1, limiting the applicability of strictlyboost invariant simulations to short times, τ (cid:46) Q − s log (1 /α s ) (19, 20, 21). In spite ofthis shortcoming, strictly boost invariant simulations of classical field dynamics form thebasis of phenomenological studies of particle production and early time dynamics in theIP-Glasma model (22, 23).During the classical evolution the field strength decreases due to the longitudinal ex-pansion, and eventually the equations of motion linearize. For times long enough τ Q s (cid:29) α − s ( Q s ). In this regime, either kinetic theory orclassical field theory can be used to simulate the evolution of the system (24, 25, 26). Inparticular, it is sensible to talk about the gluon phase space distribution, as opposed to theclassical field configuration. The initial phase space distribution of gluons f ( τ, x , p ) canbe determined from the classical simulations by evaluating the Wigner transform of equaltime two point functions of gauge fields, after fixing a physical gauge such as the CoulombGauge (see for instance ref. (27, 28, 29)). Due to the longitudinal expansion of the system, •• The First fm/c of Heavy-Ion Collisions 5 he initial phase-space distribution of the system is strongly squeezed with (cid:10) p ⊥ (cid:11) ∼ Q s and (cid:10) ( p z ) (cid:11) (cid:28) (cid:10) p ⊥ (cid:11) . This highly anisotropic initial state provides the starting point for the bottom-up scenario,which is illustrated in Fig. 2. During the first classical phase of bottom up the phase spacedistribution becomes increasingly anisotropic as time progresses.In the original bottom-up proposal, the longitudinal width of the phase space distribu-tion (cid:10) p z (cid:11) is determined by momentum diffusion, i.e. small angle scatterings amongst thehard particles. The diffusion process tries to increase the longitudinal width, but competeswith the expansion of the system. This competition leads to a scaling solution for the phasespace distribution f ( τ, p z , p ⊥ ) at late times Q s τ (cid:29)
1, where the transverse and longitudinalmomenta are of order (cid:10) p T (cid:11) ∼ Q s , (4a) (cid:10) p z (cid:11) ∼ Q s ( Q s τ ) / . (4b)During the first stage of bottom-up, the number of hard gluons per rapidity remains constant dN/dy ∼ Q s R A /α s , and thus the density of gluons (the number per volume) decreases as n h ∼ Q s /α s τ due to the expansion of the system. Based on these estimates, the phasespace density of hard modes decreases as f h ∼ α s Q s τ ) / , (5)following a pattern which is characteristic of overoccupied initial states with f h (cid:29)
1, whichwill be discussed in greater detail in Sect. 4.1 and Sect. 4.3. Analyzing eq. (5), we see thatthe phase space density becomes of order unity at a time of order Q s τ ∼ α − / s , marking theend of the first over-occupied stage. Most importantly, from this point onward the systemcan no longer be treated as a classical field, and its subsequent evolution must be analyzedwith kinetic theory.In the second stage of bottom-up, Q s τ (cid:29) α / s , radiation from the hard modes increasesthe number of soft gluons per rapidity. Ultimately this soft bath will thermalize the hardmodes giving the bottom-up equilibration scenario its name. While the soft bath is beingpopulated, the number of hard particles per volume continues to decrease due to the longi-tudinal expansion, n h ∼ Q s /α s τ . Now, however, the longitudinal width (cid:10) p z (cid:11) of these hardmodes remains constant in time, since the increase in width from (momentum) diffusion iscompensated for by the expansion of the system (cid:10) p T (cid:11) ∼ Q s , (6a) (cid:10) p z (cid:11) ∼ α s Q s . (6b)Thus, the phase space density of hard particles in the second phase decreases as f h ∼ α / s Q s τ ) , Q s τ (cid:29) α − / s , (7)and is therefore much smaller than unity. Indeed, at the end of the second phase of bottom-up, Q s τ ∼ α − / s , the phase space density of the hard modes is parametrically small, f h ∼ α s (cid:28) k
Figure 2
Schematic overview of the bottom-up thermalization showing the evolution of the phase-spacedistribution of gluons in momentum space based on kinetic theory simulations of (30). Differentregimes correspond to evolution times τ/τ hydro ≈ . , . , α s ≈ . In the last stage of the bottom-up Q s τ (cid:29) α − / s the soft bath has equilibrated, andbegins to influence the evolution of the hard particles. In this stage there is a cascade ofenergy from the scale of Q s to the soft scale scale set by the temperature of the bath. Thephysics of this process is analogous to the stopping of “jets” with momentum of order Q s in plasma (6, 31, 32) and described further in Sect. 4.2 and Sect. 4.3.The second and third stages of the bottom-up scenario are characteristic of the ther-malization of initially under-occupied systems. We will see in Sect. 4.2 that the buildup of asoft thermal bath, and cascade of energy to the infrared are to be expected in such systems.
3. QCD Kinetics: a brief review
Having qualitatively described the bottom-up picture, we will now turn to a more quanti-tative analysis of the equilibration process of the QGP in the framework of kinetic theory.Kinetic processes in the QGP are markedly different from other many-body systems of con-densed matter physics, uniquely reflecting the non-abelian and ultra-relativistic character ofthe produced quark and gluon quasi-particles. A complete leading order description of QCDkinetics (close to equilibrium) was given in (33), and was then used to compute the transportcoefficients of the QCD plasma to leading order in the strong coupling constant (34).Here we will provide a brief review of QCD kinetics to establish notation and to collectthe principal results. If not stated otherwise we will focus on pure gauge systems, and ••
Having qualitatively described the bottom-up picture, we will now turn to a more quanti-tative analysis of the equilibration process of the QGP in the framework of kinetic theory.Kinetic processes in the QGP are markedly different from other many-body systems of con-densed matter physics, uniquely reflecting the non-abelian and ultra-relativistic character ofthe produced quark and gluon quasi-particles. A complete leading order description of QCDkinetics (close to equilibrium) was given in (33), and was then used to compute the transportcoefficients of the QCD plasma to leading order in the strong coupling constant (34).Here we will provide a brief review of QCD kinetics to establish notation and to collectthe principal results. If not stated otherwise we will focus on pure gauge systems, and •• The First fm/c of Heavy-Ion Collisions 7 efer to the literature for additional details (33, 35). Further we will, at points, have toassume that the momentum distribution is isotropic; issues which arise in the descriptionof anisotropic systems (such as plasma instabilities) will be discussed briefly in Sect. 4.3.The QCD Boltzmann equation takes the form( ∂ t + v p · ∂ x ) f ( t, x , p ) = C ↔ [ f ( p )] + C ↔ [ f ( p )] , (8)where the 2 ↔ C ↔ describes collinear radiation.We further introduce two dimensionful integrals m ≡ ν g g C A d A (cid:90) p f p p , (9)2 T ∗ m ≡ ν g g C A d A (cid:90) p f p (1 + f p ) . (10)to characterize the momentum distribution. Modes of order m are considered soft, whilemodes of order T ∗ are hard. In equilibrium T ∗ is the temperature of the medium, and m isthe asymptotic mass of the gluon dispersion curve, i.e. E p = (cid:112) p + m (cid:39) | p | + m / | p | . The 2 ↔ m and screening is important, and hard collisions, where the momentum transfer isabove a cutoff scale µ ⊥ ∼ T ∗ m and screening can be neglected: C ↔ [ f ( p )] = C diff [ f ( p )] + C ↔ [ f ( p )] . (11)Hard collisions (which are conceptually straightforward) exhibit the same parametric de-pendencies as soft interactions (see e.g. (36)) and will be ignored in the estimates below.Elastic interactions with soft momentum transfers create drag and diffusion processes inmomentum space, which may be summarized by a Fokker-Planck equation. This separa-tion into hard and soft collisions was essential to an almost complete next-to-leading-ordercomputation of the shear viscosity (35).Consider a particle of momentum p (with four velocity v µ p ≡ (1 , ˆ p )) being jostled bya soft random external field A µ ( Q ) created by all other particles. The absorption rate ofthree momentum q by the field is d Γ > el ( ˆ p ) d q = g C A (cid:90) dq π v µ p v ν p (cid:104)(cid:104) A µ ( Q ) ( A ν ( Q )) ∗ (cid:105)(cid:105) > πδ ( v p · Q ) , (12)where the δ -function stems from energy conservation, 2 πδ ( q + E p − q − E p ) (cid:39) πδ ( v p · Q ).Statistical fluctuations of the gauge field fluctuations are given by (cid:104)(cid:104) A µ ( Q )( A ν ( Q )) ∗ (cid:105)(cid:105) > = G µαR ( Q ) Π >αβ ( Q ) ( G βνR ( Q )) ∗ , (13)where the Wightman self energy reads(Π > ( Q )) αβ = ν g g C A d A (cid:90) k v α k v β k f ( k )(1 + f ( k + q )) 2 πδ ( v k · Q ) , (14) We follow standard notation, where d A = N c − C A = N c is its Casimir. By ν g = 2 d A we denote the number of gluonic degrees of freedom. Phase spaceintegrals are abbreviated as (cid:82) p ≡ (cid:82) d p/ (2 π ) , as is the phase space density f p ≡ f ( t, x , p ). ere G R ( Q ) ∼ /Q is the hard thermal loop retarded response function (37), whichcan only be worked out in closed form for isotropic systems. In the limit of small q thepopulation factors in eq. (14) become f ( k )(1 + f ( k )), and the correlator in eq. (13) has asimple interpretation – it is the correlation amongst the gauge fields A = ( G R ( Q )) · gv k produced by random fluctuations of the phase space density δf ( t, x , k ), which have theusual equal time Bose-Einstein statistics (38) (cid:104)(cid:104) δf ( t, x , k ) δf ( t, x (cid:48) , k ) (cid:105)(cid:105) = f ( t, x , k ) (1 + f ( t, x , k )) δ ( x − x (cid:48) ) (2 π ) δ ( k − k (cid:48) ) . (15)The absorption rate in eq. (12) gives the rate that momentum q is taken from theparticle and given to the bath. Similarly, the emission rate takes the same form as eq. (12)but replaces the self energy Π > with(Π < ( Q )) αβ = ν g g C A d A (cid:90) k v α k v β k f ( k + q )(1 + f ( k )) 2 πδ ( v k · Q ) , (16)such that at small q the emission and absorption rates are equal, and it is the symmetriccorrelator Γ el =(Γ > el + Γ < el ) / > ( Q ) − Π < ( Q )) αβ = ν g g C A d A (cid:90) k v α k v β k q i ∂f ( k ) ∂k i πδ ( v k · Q ) , (17)=2 q m (cid:90) d Ω4 π v α k v β k πδ ( v k · Q ) , (18)where in passing to the last line we have assumed that the system is isotropic, ∂f/∂k i = f (cid:48) ( k )ˆ k i , allowing us to perform an integration by parts.The evolution of the system due to soft scattering is a competition between the emissionand absorption rates ∂ t f p + v p · ∂ x f p = (cid:90) d q (cid:18) dΓ < el ( ˆ p )d q f p − q (1 + f p ) − dΓ > el ( ˆ p )d q f p (1 + f p − q ) (cid:19) . (19)We will now generally assume that the distribution is isotropic which simplifies the analysisof momentum diffusion. Expanding in powers of the momentum transfer q (which is smallcompared to the momentum p of the hard particle), we see that the contribution of smallangle elastic processes to the Boltzmann equation (8) takes the form of a Fokker-Planckequation C diff [ f ( p )] = η i ( ˆ p ) ∂∂p i ( f p (1 + f p )) + ˆ q ij ( ˆ p ) ∂ f p ∂p i ∂p j , (20)where the drag and diffusion coefficients are given by η i = (cid:90) d q (cid:18) dΓ > el ( ˆ p )d q − dΓ < el ( ˆ p )d q (cid:19) q i , (21)ˆ q ij ( ˆ p ) = (cid:90) d q (cid:18) dΓ el ( ˆ p )d q (cid:19) q i q j . (22)Specifically for isotropic systems these coefficients can be decomposed as η i ( ˆ p ) = η ˆ p i , q ij ( ˆ p ) = ˆ q L ˆ p i ˆ p j + ˆ q (cid:16) δ ij − ˆ p i ˆ p j (cid:17) , (23) ••
Having qualitatively described the bottom-up picture, we will now turn to a more quanti-tative analysis of the equilibration process of the QGP in the framework of kinetic theory.Kinetic processes in the QGP are markedly different from other many-body systems of con-densed matter physics, uniquely reflecting the non-abelian and ultra-relativistic character ofthe produced quark and gluon quasi-particles. A complete leading order description of QCDkinetics (close to equilibrium) was given in (33), and was then used to compute the transportcoefficients of the QCD plasma to leading order in the strong coupling constant (34).Here we will provide a brief review of QCD kinetics to establish notation and to collectthe principal results. If not stated otherwise we will focus on pure gauge systems, and •• The First fm/c of Heavy-Ion Collisions 7 efer to the literature for additional details (33, 35). Further we will, at points, have toassume that the momentum distribution is isotropic; issues which arise in the descriptionof anisotropic systems (such as plasma instabilities) will be discussed briefly in Sect. 4.3.The QCD Boltzmann equation takes the form( ∂ t + v p · ∂ x ) f ( t, x , p ) = C ↔ [ f ( p )] + C ↔ [ f ( p )] , (8)where the 2 ↔ C ↔ describes collinear radiation.We further introduce two dimensionful integrals m ≡ ν g g C A d A (cid:90) p f p p , (9)2 T ∗ m ≡ ν g g C A d A (cid:90) p f p (1 + f p ) . (10)to characterize the momentum distribution. Modes of order m are considered soft, whilemodes of order T ∗ are hard. In equilibrium T ∗ is the temperature of the medium, and m isthe asymptotic mass of the gluon dispersion curve, i.e. E p = (cid:112) p + m (cid:39) | p | + m / | p | . The 2 ↔ m and screening is important, and hard collisions, where the momentum transfer isabove a cutoff scale µ ⊥ ∼ T ∗ m and screening can be neglected: C ↔ [ f ( p )] = C diff [ f ( p )] + C ↔ [ f ( p )] . (11)Hard collisions (which are conceptually straightforward) exhibit the same parametric de-pendencies as soft interactions (see e.g. (36)) and will be ignored in the estimates below.Elastic interactions with soft momentum transfers create drag and diffusion processes inmomentum space, which may be summarized by a Fokker-Planck equation. This separa-tion into hard and soft collisions was essential to an almost complete next-to-leading-ordercomputation of the shear viscosity (35).Consider a particle of momentum p (with four velocity v µ p ≡ (1 , ˆ p )) being jostled bya soft random external field A µ ( Q ) created by all other particles. The absorption rate ofthree momentum q by the field is d Γ > el ( ˆ p ) d q = g C A (cid:90) dq π v µ p v ν p (cid:104)(cid:104) A µ ( Q ) ( A ν ( Q )) ∗ (cid:105)(cid:105) > πδ ( v p · Q ) , (12)where the δ -function stems from energy conservation, 2 πδ ( q + E p − q − E p ) (cid:39) πδ ( v p · Q ).Statistical fluctuations of the gauge field fluctuations are given by (cid:104)(cid:104) A µ ( Q )( A ν ( Q )) ∗ (cid:105)(cid:105) > = G µαR ( Q ) Π >αβ ( Q ) ( G βνR ( Q )) ∗ , (13)where the Wightman self energy reads(Π > ( Q )) αβ = ν g g C A d A (cid:90) k v α k v β k f ( k )(1 + f ( k + q )) 2 πδ ( v k · Q ) , (14) We follow standard notation, where d A = N c − C A = N c is its Casimir. By ν g = 2 d A we denote the number of gluonic degrees of freedom. Phase spaceintegrals are abbreviated as (cid:82) p ≡ (cid:82) d p/ (2 π ) , as is the phase space density f p ≡ f ( t, x , p ). ere G R ( Q ) ∼ /Q is the hard thermal loop retarded response function (37), whichcan only be worked out in closed form for isotropic systems. In the limit of small q thepopulation factors in eq. (14) become f ( k )(1 + f ( k )), and the correlator in eq. (13) has asimple interpretation – it is the correlation amongst the gauge fields A = ( G R ( Q )) · gv k produced by random fluctuations of the phase space density δf ( t, x , k ), which have theusual equal time Bose-Einstein statistics (38) (cid:104)(cid:104) δf ( t, x , k ) δf ( t, x (cid:48) , k ) (cid:105)(cid:105) = f ( t, x , k ) (1 + f ( t, x , k )) δ ( x − x (cid:48) ) (2 π ) δ ( k − k (cid:48) ) . (15)The absorption rate in eq. (12) gives the rate that momentum q is taken from theparticle and given to the bath. Similarly, the emission rate takes the same form as eq. (12)but replaces the self energy Π > with(Π < ( Q )) αβ = ν g g C A d A (cid:90) k v α k v β k f ( k + q )(1 + f ( k )) 2 πδ ( v k · Q ) , (16)such that at small q the emission and absorption rates are equal, and it is the symmetriccorrelator Γ el =(Γ > el + Γ < el ) / > ( Q ) − Π < ( Q )) αβ = ν g g C A d A (cid:90) k v α k v β k q i ∂f ( k ) ∂k i πδ ( v k · Q ) , (17)=2 q m (cid:90) d Ω4 π v α k v β k πδ ( v k · Q ) , (18)where in passing to the last line we have assumed that the system is isotropic, ∂f/∂k i = f (cid:48) ( k )ˆ k i , allowing us to perform an integration by parts.The evolution of the system due to soft scattering is a competition between the emissionand absorption rates ∂ t f p + v p · ∂ x f p = (cid:90) d q (cid:18) dΓ < el ( ˆ p )d q f p − q (1 + f p ) − dΓ > el ( ˆ p )d q f p (1 + f p − q ) (cid:19) . (19)We will now generally assume that the distribution is isotropic which simplifies the analysisof momentum diffusion. Expanding in powers of the momentum transfer q (which is smallcompared to the momentum p of the hard particle), we see that the contribution of smallangle elastic processes to the Boltzmann equation (8) takes the form of a Fokker-Planckequation C diff [ f ( p )] = η i ( ˆ p ) ∂∂p i ( f p (1 + f p )) + ˆ q ij ( ˆ p ) ∂ f p ∂p i ∂p j , (20)where the drag and diffusion coefficients are given by η i = (cid:90) d q (cid:18) dΓ > el ( ˆ p )d q − dΓ < el ( ˆ p )d q (cid:19) q i , (21)ˆ q ij ( ˆ p ) = (cid:90) d q (cid:18) dΓ el ( ˆ p )d q (cid:19) q i q j . (22)Specifically for isotropic systems these coefficients can be decomposed as η i ( ˆ p ) = η ˆ p i , q ij ( ˆ p ) = ˆ q L ˆ p i ˆ p j + ˆ q (cid:16) δ ij − ˆ p i ˆ p j (cid:17) , (23) •• The First fm/c of Heavy-Ion Collisions 9 nd the scalar coefficients η, ˆ q L , ˆ q can be evaluated as (see (39) for a review), η = g C A m π log (cid:18) µ ⊥ m (cid:19) , (24a)ˆ q L = g C A (2 T ∗ m )8 π log (cid:18) µ ⊥ m (cid:19) , (24b)ˆ q = g C A (2 T ∗ m )4 π log (cid:18) µ ⊥ m (cid:19) . (24c)Similarly, the elastic scattering rate for kicks transverse to the direction of the particle canalso be evaluated in closed form yielding(2 π ) d Γ el d q ⊥ = g C A T ∗ (cid:18) q ⊥ − q ⊥ + 2 m (cid:19) . (25)Although the Fokker-Planck coefficients in eq. (24a) depend on the cutoff scale µ ⊥ , thetime the evolution of the system is independent of µ ⊥ , when both the hard collisions andthe Fokker-Planck evolution are taken into account (40). We finally note that from eq. (25)and eq. (24c), the elastic scattering rate is of orderΓ el ∼ (cid:90) ∼ m d q ⊥ dΓ el d q ⊥ ∼ ˆ qm , (26)which will be used repeatedly when estimating the rate of collinear radiation described inthe next section. Elastic scatterings of ultra-relativistic particles induce collinear radiation as the chargedparticles are accelerated by the random kicks from the plasma. A massless gluon withmomentum P = p + k can split into two particles with momentum fractions z and ¯ z ≡ (1 − z ),where p = z P and k = ¯ z P respectively. These radiative process should be incorporatedinto the Boltzmann equation at leading order (6, 33). Denoting the rate for this process asdΓ inel ( P ) / d z , the contribution to the Boltzmann equation can be written as C ↔ [ f ( p )] = ν g (cid:90) P (cid:90) dz dΓ inel ( P )d z (2 π ) ν g δ (3) ( p − z P ) × [ f ( P )(1 + f ( z p ))(1 + f (¯ z P )) − f ( z P ) f (¯ z P )(1 + f ( P ))] − (cid:90) dz d Γ inel ( p ) dz × [ f ( p )(1 + f ( z p ))(1 + f (¯ z p )) − f ( z p ) f (¯ z p )(1 + f ( p ))] , (28)and we will now briefly describe the characteristic features of the splitting rate. Our notation for inelastic splitting rate follows (41, 42). Arnold, Moore, and Yaffe use adifferent symbol γ ggg ( p (cid:48) , p , k ) (33), which is related to the rate used here through d Γ inel ( P ) dz = (2 π ) ν g | P | γ ggg ( P , z P , (1 − z ) P ) . (27)
10 S. Schlichting and D. Teaney n the splitting process the energy difference between the incoming and outgoing statesis δE = E p + E k − E p + k (cid:39) h P z (1 − z ) + m P z + m P (1 − z ) − m P , (29)where h ≡ z k ⊥ − (1 − z ) p ⊥ is essentially the transverse momentum of the softest fragment. Inwriting eq. (29) we have expanded the quasiparticle energy for small transverse momentum, E p (cid:39) p z +( m + p ⊥ ) / p . Since the Hamiltonian time evolution of the system involves phasesof the form e − iδEt , the splitting process is only completed on a time scale t form ≡ δE . (30)which defines an important timescale for collinear radiation, namely the formation time .For highly energetic particles the formation time can become long compared to the timebetween elastic collisions. In this regime multiple scattering will suppress the emission ofgluon radition, and this suppression is known as the Landau-Pomenanchuk-Migdal (LPM)effect.Let us estimate the energy ω LPM when the LPM effect becomes operative, i.e. when t form Γ el ∼
1. To this end, consider a splitting process with z (cid:28) h (cid:39) p ⊥ and p = zP ∼ ω LPM . In this regime the formation time is of order t form ∼ pp ⊥ ∼ ω LPM m , (31)were we have estimated p ⊥ ∼ m as the typical momentum associated with a single elasticscattering event. Since the elastic scattering rate is of order Γ el ∼ ˆ q/m , we find ω LPM ∼ m ˆ q . (32)For high energy particles the formation time becomes much longer than Γ − . In this limitthe accumulated transverse momentumg grows as h ∼ ˆ q t form (cid:29) m , and thus using eq. (30)and eq. (29) we find the following estimate for the formation time t form ∼ (cid:115) Pz (1 − z )ˆ q . (33)For ω (cid:38) ω LPM the radiation rate must account for the multiple scatterings that happenduring the formation time of the radiation. Conversely, in the Bethe-Heitler (BH) limit ω (cid:28) ω LPM , the interference between the scattering events can be neglected, and each scat-tering has a probability of order α to radiate a gluon with momentum fraction z disributedaccording to the splitting function P g → g ( z ). Since the scattering rate is Γ el ∼ ˆ q/m , thetotal splitting rate in the BH limit is of orderdΓ BH inel ( p )d z ∼ α P soft g → g ( z ) ˆ qm . (34) Generally the splitting function for g ↔ gg is given by P g → g ( z ) = C A z +(1 − z ) z (1 − z ) . Howeverwe will frequently approximate P g → g ( z ) by its soft limit P soft g → g ( z ) = C A z (1 − z ) . ••
1. To this end, consider a splitting process with z (cid:28) h (cid:39) p ⊥ and p = zP ∼ ω LPM . In this regime the formation time is of order t form ∼ pp ⊥ ∼ ω LPM m , (31)were we have estimated p ⊥ ∼ m as the typical momentum associated with a single elasticscattering event. Since the elastic scattering rate is of order Γ el ∼ ˆ q/m , we find ω LPM ∼ m ˆ q . (32)For high energy particles the formation time becomes much longer than Γ − . In this limitthe accumulated transverse momentumg grows as h ∼ ˆ q t form (cid:29) m , and thus using eq. (30)and eq. (29) we find the following estimate for the formation time t form ∼ (cid:115) Pz (1 − z )ˆ q . (33)For ω (cid:38) ω LPM the radiation rate must account for the multiple scatterings that happenduring the formation time of the radiation. Conversely, in the Bethe-Heitler (BH) limit ω (cid:28) ω LPM , the interference between the scattering events can be neglected, and each scat-tering has a probability of order α to radiate a gluon with momentum fraction z disributedaccording to the splitting function P g → g ( z ). Since the scattering rate is Γ el ∼ ˆ q/m , thetotal splitting rate in the BH limit is of orderdΓ BH inel ( p )d z ∼ α P soft g → g ( z ) ˆ qm . (34) Generally the splitting function for g ↔ gg is given by P g → g ( z ) = C A z +(1 − z ) z (1 − z ) . Howeverwe will frequently approximate P g → g ( z ) by its soft limit P soft g → g ( z ) = C A z (1 − z ) . •• The First fm/c of Heavy-Ion Collisions 11 ore generally emissions radiated within a formation time will destructively interfere, andthe net emission rate is determined by solving an integral equation. This rate takes theform (43, 44, 6, 33, 41) d Γ inel ( P )d z = α s P g → g ( z ) (cid:90) d h (2 π ) h · Re f ( h )(2 P z (1 − z )) , (35)where the integral in this equation has units (time) − . The function f ( h ) (which encodesthe current-current statistical correlation function) satisfies an integral equation of the form2 h = i δE ( h ) f ( h ) + (cid:90) d q ⊥ d Γ el d q ⊥ (cid:8) [ f ( h ) − f ( h + q ⊥ )]+ [ f ( h ) − f ( h + z q ⊥ )] + [ f ( h ) − f ( h + (1 − z ) q ⊥ )] (cid:9) . (36)To analyze this equation, let us take the Bethe-Heitler limit when the radiation is soft, z (cid:28) ω (cid:28) ω LPM , so that the formation time is small compared to the elastic scatteringrate, δE (cid:29) Γ el . In this regime we can solve eq. (36) by iteration, f = f (0) + f (1) + . . . , with f (0) ( h ) = − i h /δE ( h ). Physically this expansion corresponds to the number of collisions,with f (1) determining the emission rate from one collision and so on. After straightforwardalgebra one finds dΓ BH inel ( p )d z = 2 α s P soft g → g ( z ) (cid:90) d p ⊥ (2 π ) (cid:90) d q ⊥ d Γ el d q ⊥ (cid:18) p ⊥ p ⊥ + m − p ⊥ + q ⊥ ( p ⊥ + q ⊥ ) + m (cid:19) . (37)The large p ⊥ limit of this rate is known as the Gunion-Bertsch formula (45)(2 π ) dΓ BH inel ( P )d z d p ⊥ = 2 α s P soft g → g ( z ) ˆ qp ⊥ . (38)To estimate the total rate one can integrate this expression over p ⊥ down to a scale p ⊥ ∼ m yielding the Bethe-Heitler estimate given earlier in eq. (34) .In the opposite limit ω (cid:29) ω LPM we can also find an approximate solution to eq. (36)known as the harmonic oscillator approximation. Since for ω (cid:29) ω LPM the transversemomentum h acquired over the formation time is large compared to the typical momentumtransfer q ⊥ aquired in a single scattering q ⊥ ∼ m , one can expand the differences f ( h ) − f ( h + q ⊥ ) for small q ⊥ , which transforms (36) into a partial differential equation2 h = iδE ( h ) f ( h ) − z + (1 − z ) q δ ij ⊥ ∂ ∂h i ∂h j f ( h ) (39)By approximating δE ( h ) (cid:39) h Pz (1 − z ) and Fourier transforming with respect to h (with b conjugate to h ), one obtains a Schr¨odinger-like equation for a particle with an effectivemass M = P z (1 − z ) in an imaginary harmonic potential V ( b ) = − i Mω b with oscillationfrequency ω = ˆ q z +(1 − z ) z (1 − z ) P . Solving this equation, one finds that the final emission rate Note that we have somewhat cavalierly shifted the integration variable p ⊥ → p ⊥ + q ⊥ to re-write p ⊥ δE ( p ⊥ ) → (cid:18) p ⊥ δE ( p ⊥ ) + ( p ⊥ + q ⊥ ) δE ( p ⊥ + q ⊥ ) (cid:19) in order to write the integrand as a perfect square,which naturally appears in diagrammatic calculations of the single scattering rates (40).
12 S. Schlichting and D. Teaney s proportional to ω ∼ t − , which in the soft limit ( z (cid:28)
1) yields d Γ LPM inel ( P ) dz = α s π P soft g → g ( z ) (cid:115) ˆ qP z (1 − z ) . (40)General expressions for the emission rates involving multiple species are given in (46, 42) inthe same notation used here. Comparing eq. (40) with the Bethe-Heitler limit of eq. (34),shows that the emission rate is controlled by the inverse of the formation time 1 /t form ratherthan the elastic scattering rate ∼ ˆ q/m in eq. (34), suppressing the emission of radiation athigh energies.
4. Basics of weak coupling thermalization
Now that we have outlined the basic physics of QCD kinetics, we will illustrate key featuresof the equilibration process in homogenous isotropic systems where a detailed understand-ing of the dynamics has been gained in a series of studies (36, 47, 48, 49). Since theequilibration dynamics crucially depends on the properties of the initial state, it useful todistinguish between systems which are initially far from equilibrium, and systems whichare initially close equilibrium. While in the latter case, one expects a direct relaxation ofthe system to equilibrium governed by an equilibrium rate, the situation is more compli-cated for systems which are initially far from equilibrium, and various kinds of phenomenacan occur en-route towards thermal equilibrium. Nevertheless a general characterization ofthe equilibration process can be achieved for broad classes of far from equilibrium initialconditions. Specifically, for homogenous and isotropic systems one needs to distinguishbetween overoccupied systems, i.e. systems in which the energy is initially carried by alarge number of low energy degrees of freedom, and underoccupied systems, i.e. systemsin which the energy is carried by a small number of very high energy degrees of freedom.As we have emphasized in the previous section the first stage of the bottom-up scenariocorresponds to the “over-occupied” case, while the second and third stages correspond tothe under-occuppied case.
We first consider a system where the initial energy density is carried by a large number oflow energy degrees of freedom, i.e. if the quasi-particle energy is E p ∼ Q , then the energydensity is e ∼ f Q where f (cid:29) e eq ∼ T is carried bya smaller number of modes with f ∼ E p ∼ T . Since energy is conserved during theevolution, the final temperatue T ∼ Q f / at the end of the equilibration process is muchlarger than Q . Because of this large scale separation between Q and T , the redistribution ofenergy from low energy modes to high energy modes is then a classic problem of turbulenceknown as a direct energy cascade discussed in the next section (50). The initial evolution of overoccu-pied plasmas can be equivalently described in terms of classical fields or weakly interactingquasi particles, due to an overlap in their respective range of validity (24, 25, 26). For thisreason the initial evolution can either be studied using classical-statistical simulations of ••
We first consider a system where the initial energy density is carried by a large number oflow energy degrees of freedom, i.e. if the quasi-particle energy is E p ∼ Q , then the energydensity is e ∼ f Q where f (cid:29) e eq ∼ T is carried bya smaller number of modes with f ∼ E p ∼ T . Since energy is conserved during theevolution, the final temperatue T ∼ Q f / at the end of the equilibration process is muchlarger than Q . Because of this large scale separation between Q and T , the redistribution ofenergy from low energy modes to high energy modes is then a classic problem of turbulenceknown as a direct energy cascade discussed in the next section (50). The initial evolution of overoccu-pied plasmas can be equivalently described in terms of classical fields or weakly interactingquasi particles, due to an overlap in their respective range of validity (24, 25, 26). For thisreason the initial evolution can either be studied using classical-statistical simulations of •• The First fm/c of Heavy-Ion Collisions 13 he non-linear gauge field dynamics (see e.g. (48)), or using the numerical simulations andanalytic considerations of kinetic theory (36, 47).It was found that the initial evolution of overoccupied systems proceeds via a quasi-stationary state referred to as a non-thermal fixed point (NTFP). Here the dynamics be-comes insensitive to the details of the initial conditions after a short time, and the evolutionfollows a self-similar scaling behavior (51, 52, 48). Indeed, the phase-space density f ( t, p )in this regime evolves with the scaling form f ( t, p ) = ( Qt ) α f S (cid:18) ( Qt ) β pQ (cid:19) , (41)which is characteristic for non-stationary turbulent processes (50) and the scaling form ineq. (41) describes a direct energy cascade , i.e. the transport of energy from low momentumto high momentum excitations necessary to achieve thermalization.The scaling exponents α, β (which will be negative) determine the increase of the char-acteristic momentum scale p max ( t ) ∼ Q ( Qt ) − β , and the simultaneous decrease of the occu-pancy of hard excitations f ( t, p ∼ p max ( t )) ∼ ( Qt ) α (see Fig. 3). These scaling exponentscan be determined from a straightforward scaling analysis of the underlying kinetic equa-tions (47, 36, 51, 48) following well established techniques in the context of weak waveturbulence (50). One immediate constraint on the scaling exponents α, β comes from therequirement of energy conservation e ( t ) = (cid:90) p E p f ( t, p ) = const , (42)which for a self-similar evolution of the form in eq. (41) gives rise to a scaling relation α − β = 0 . (43)A second scaling relation can be inferred from a scaling analysis of the kinetic equation.Even though the full scaling analysis of all leading order kinetic processes is somewhatcomplicated, the essence can be understood by considering as an example small angle elasticprocesses, whose contribution to the collision integral in eq. (20) is of the form of a Fokker-Planck equation, where the drag coefficient η ( t ) and momentum diffusion coefficient ˆ q ( t )are of order (see eq. (24) and eq. (9)) η ( t ) ∼ α s (cid:90) p f ( t, p ) p , (44a)ˆ q ( t ) ∼ α s (cid:90) p f ( t, p ) (1 + f ( t, p )) . (44b)With the scaling ansatz of eq. (41) in the high-occupancy regime f ( t, p ) (cid:29)
1, these quan-tities scale (up to logarithmic corrections) as η ( t ) ∼ ( Qt ) α − β α s Q (cid:90) q f S ( q ) q , (45)ˆ q ( t ) ∼ ( Qt ) α − β α s Q (cid:90) q f S ( q ) , (46)under the self-similar evolution of the system. Based on this analysis one can establish ascaling behavior of the collision integral C diff [ f ( t, p )] = ( Qt ) α − β C diff [ f S ( Q )] , (47)
14 S. Schlichting and D. Teaney s e l f - s i m i l a r d i r ec t e n e r g y c a s c ad e f Evolution of over-occupied plasma T1f(t,p) p max (t)=Q(Qt) - β Q i n ve r s e e n e r g y c a s c ad e f Evolution of under-occupied plasmaT1 p soft (t) s o f t t h e r m a l b a t h p split (t)f(t,p) Figure 3
Illustration of the thermalization process in over-occupied and under-occupied systems,summarizing the results of classical-statistical field simulations(52, 51, 48) and kinetic theorysimulations (53, 49). which also extends to large angle elastic and inelastic processes (47, 36, 51, 48). By matchingthe time dependence on the r.h.s of Eq. (47) with that of the l.h.s. of the Boltzmannequation, one infers the dynamical scaling relation α − α − β , (48)which along with Eq. (43) uniquely determines the exponents. Strikingly, the scaling anal-ysis of the kinetic equations also reveals the universal nature of the dynamical scalingexponents, which are insensitive to microscopic details of the underlying theory and takethe values α = − / β = − / SU ( N c ) gauge theories in d = 3 dimensions (47, 36).These are in line with classical-statistical field simulations of SU (2) and SU (3) Yang-Millsplasmas (51, 48, 54). ••
Illustration of the thermalization process in over-occupied and under-occupied systems,summarizing the results of classical-statistical field simulations(52, 51, 48) and kinetic theorysimulations (53, 49). which also extends to large angle elastic and inelastic processes (47, 36, 51, 48). By matchingthe time dependence on the r.h.s of Eq. (47) with that of the l.h.s. of the Boltzmannequation, one infers the dynamical scaling relation α − α − β , (48)which along with Eq. (43) uniquely determines the exponents. Strikingly, the scaling anal-ysis of the kinetic equations also reveals the universal nature of the dynamical scalingexponents, which are insensitive to microscopic details of the underlying theory and takethe values α = − / β = − / SU ( N c ) gauge theories in d = 3 dimensions (47, 36).These are in line with classical-statistical field simulations of SU (2) and SU (3) Yang-Millsplasmas (51, 48, 54). •• The First fm/c of Heavy-Ion Collisions 15 eyond the dynamics of energy transport, various perturbative and non-perturbativeproperties of the NTFPs of SU ( N c ) Yang-Mills theory have been investigated based onclassical-statistical lattice simulations (55, 54, 56) and show how the electric and magneticsectors of kinetic theory emerge at late time. Eventually, the self-similar evolution breaks down when the energyhas been transferred from the initial momentum scale p max ( t = 0) ∼ Q all the way to theequilibrium temperature p max ( t eq ) ∼ T (47, 36). Using the scaling exponent β and theinitial occupancy f ∼ /α s , the self-similar cascade ends when t ∼ t eq ∼ α − s f − / Q − ∼ α − s T − . (49)At the end of the cascade, the phase-space occupancies of hard modes f ( t, p max ( t )) alsobecomes of order unity, and the system is no longer parametrically far from equilibrium.The relevant scattering rates decrease over the course of the cascade, Γ( t ) ∼ ˆ q ( t ) p − ( t ) ∼ α s Q ( Qt ) − , and the final approach to equilibrium is ultimately controlled by an equilibriumtransport time scale, ∼ α − s T − . This time scale is parametrically of the same order asthe time scale for the turbulent transport of energy given in eq. (49). While the finalapproach to equilibrium is outside the range of validity of classical-statistical simulations,it can be investigated further based on numerical simulations in kinetic theory (53, 49),which provide concrete, rather than parametric, estimates of the thermalization time, t eq ≈ . α − s N − c T − (49). We now consider the opposite case where the initial energy density e ∼ Q is carried bya small number f (cid:28) E p ∼ Q , and note thatthis setup is reminiscent of a high-energy jet carrying a significant fraction of the energyof the system. While the final equilibrium temperature can again be inferred using energyconservation as T ∼ f / Q , the hierarchy of scales is now inverted with T (cid:28) Q . Sincethe equilibrium temperature T is much smaller than the characteristic momentum scale Q ,the thermalization process now requires a re-distribution of energy from high energy to lowenergy degrees of freedom.Eventually the re-distribution of energy is achieved by an inverse energy cascade throughmultiple radiative branchings of the high energy particles (6, 36, 32). However, before theinverse cascade can be established, a small fraction of the energy must be transferred tolow energy modes by direct emission of soft radiation. As discussed in Sect. 4.2.1, these lowenergy modes thermalize quickly, creating of a soft thermal bath and setting the stage forthe inverse energy cascade described in Sect. 4.2.2. Let us analyze how the soft bathis created. Following (36) there is a competition between direct radiation from the hardmodes, which tends to populate the soft bath, and momentum diffusion which tends topush the typical momentum scale of the bath to higher momentum. As we will show below,direct radiation initially dominates and over populates the bath. Then, as the LPM effectsets in and suppresses additional radiation, the soft bath reaches an occupancy of orderunity with an equilibrium temperature T soft ( t ).Initially, elastic scattering processes amongst the hard modes occur relatively frequently,
16 S. Schlichting and D. Teaney ith a rate of order Γ el ∼ ˆ qm ∼ α s Q , where we have estimated ˆ q and m from the distribution of hard particlesˆ q hard ∼ α s (cid:90) p f p (1 + f p ) ∼ α s f Q , (50a) m ∼ α s (cid:90) p f p p ∼ α s f Q . (50b)These elastic scatterings induce soft and collinear radiation, and it is these processes whichare responsible for creating the soft bath. From the first line of eq. (28), the rate at whichsoft particles with momentum p are produced by the hard particles with momentum P ∼ Q is initially ∂f ( t, p ) ∂t (cid:39) ν g (cid:90) P (cid:90) dz dΓ inel ( P )d z (2 π ) ν g δ (3) ( p − z P ) f ( P ) (cid:16) f ( P ) (cid:17) . (51)Note that this rate is independent of the soft phase space density f ( t, p ) due to a cancellationbetween the gain and loss terms (36).The radiated soft fragments are of course more susceptible to elastic scattering processes,and have the chance to equilibrate via both elastic scatterings and inelastic processes, givingrise to a dynamical scale p soft ( t ) ∼ (cid:112) ˆ q ( t ) t ∼ α s f / Q ( Qt ) / . (52)Soft fragments below p soft ( t ) have an effective temperature T ∗ soft ( t ) (defined precisely below)characterizing the occupancy of these modes.As we will now estimate, the phase space densities become initially overoccuppied as thesoft bath is built up. This happens because at early times the particles are copiously pro-duced via Bethe-Heitler radiation, and do not have time to increase p soft through diffusion.The Bethe-Heitler approximation is appropriate here because p soft (cid:28) ω LPM as discussed inSect. 3.2. The occupancy of the soft sector can be estimated from the amount of energy e soft radiated into this sector and p soft ( t ). The radiated energy is of order e soft ( t ) ∼ (cid:90) t dt (cid:90) p max ( t ) p E p ∂f ( t, | p | ) ∂t ∼ e hard (cid:90) p soft ( t ) /Q dz z dΓ inel ( Q )d z t , (53)which, with the Bethe-Heitler estimate for dΓ BH inel / d z from eq. (34), yields e soft ( t ) ∼ α s e hard (ˆ q ( t ) t ) m p soft ( t ) Q . (54)Using the estimates for ˆ q and m in eq. (50), one finds that the effective temperature of thesoft sector is given by T ∗ soft ( t ) ≡ e soft ( t )( p soft ( t )) ∼ Q , (55)and thus, since T ∗ soft ( t ) is much larger than the characteristic momentum scale p soft ( t ), thesystem is initially over occupied for a short period of time( Qt ) / (cid:46) α − s f − / . (56) ••
16 S. Schlichting and D. Teaney ith a rate of order Γ el ∼ ˆ qm ∼ α s Q , where we have estimated ˆ q and m from the distribution of hard particlesˆ q hard ∼ α s (cid:90) p f p (1 + f p ) ∼ α s f Q , (50a) m ∼ α s (cid:90) p f p p ∼ α s f Q . (50b)These elastic scatterings induce soft and collinear radiation, and it is these processes whichare responsible for creating the soft bath. From the first line of eq. (28), the rate at whichsoft particles with momentum p are produced by the hard particles with momentum P ∼ Q is initially ∂f ( t, p ) ∂t (cid:39) ν g (cid:90) P (cid:90) dz dΓ inel ( P )d z (2 π ) ν g δ (3) ( p − z P ) f ( P ) (cid:16) f ( P ) (cid:17) . (51)Note that this rate is independent of the soft phase space density f ( t, p ) due to a cancellationbetween the gain and loss terms (36).The radiated soft fragments are of course more susceptible to elastic scattering processes,and have the chance to equilibrate via both elastic scatterings and inelastic processes, givingrise to a dynamical scale p soft ( t ) ∼ (cid:112) ˆ q ( t ) t ∼ α s f / Q ( Qt ) / . (52)Soft fragments below p soft ( t ) have an effective temperature T ∗ soft ( t ) (defined precisely below)characterizing the occupancy of these modes.As we will now estimate, the phase space densities become initially overoccuppied as thesoft bath is built up. This happens because at early times the particles are copiously pro-duced via Bethe-Heitler radiation, and do not have time to increase p soft through diffusion.The Bethe-Heitler approximation is appropriate here because p soft (cid:28) ω LPM as discussed inSect. 3.2. The occupancy of the soft sector can be estimated from the amount of energy e soft radiated into this sector and p soft ( t ). The radiated energy is of order e soft ( t ) ∼ (cid:90) t dt (cid:90) p max ( t ) p E p ∂f ( t, | p | ) ∂t ∼ e hard (cid:90) p soft ( t ) /Q dz z dΓ inel ( Q )d z t , (53)which, with the Bethe-Heitler estimate for dΓ BH inel / d z from eq. (34), yields e soft ( t ) ∼ α s e hard (ˆ q ( t ) t ) m p soft ( t ) Q . (54)Using the estimates for ˆ q and m in eq. (50), one finds that the effective temperature of thesoft sector is given by T ∗ soft ( t ) ≡ e soft ( t )( p soft ( t )) ∼ Q , (55)and thus, since T ∗ soft ( t ) is much larger than the characteristic momentum scale p soft ( t ), thesystem is initially over occupied for a short period of time( Qt ) / (cid:46) α − s f − / . (56) •• The First fm/c of Heavy-Ion Collisions 17 he radiated soft excitations will ultimately contribute to screening and scattering pro-cesses. While at early times these contributions are negligible, their contributions increaseas a function of time according to m ( t ) ∼ α s T ∗ soft ( t ) p soft ( t ) ∼ m ( Qt ) / f − / α − s , (57a)ˆ q soft ( t ) ∼ α s ( T ∗ soft ( t )) p soft ( t ) ∼ ˆ q hard ( Qt ) / f − / α − s , (57b)and thus for ( Qt ) / (cid:38) α − s f − / they become of the same order as the contributions fromthe hard sector, and the systems enters the second stage of the thermalization process.For ( Qt ) / (cid:38) α − s f − / , the radiative dynamics continues in a similar fashion, but nowthe soft and hard sectors now give comparable contributions to elastic scattering, while thescreening is dominated by the soft sector. The emission of soft radiation at the characteristicscale p soft ( t ) now suffers from LPM suppression as now p soft ( t ) has become of order of ω LPM .Substituting eq. (40) in eq. (53), the amount of energy radiated directly into soft modes p ∼ p soft ( t ) is now given by e soft ( t ) ∼ α s e hard (cid:115) ˆ q ( t ) t Q (cid:115) p soft ( t ) Q , (58)which along with the consistency relations e soft ∼ T ∗ soft ( t ) p ( t ) , p soft ( t ) ∼ (cid:112) ˆ q ( t ) t , ˆ q ( t ) ∼ ˆ q soft ( t ) ∼ α s T ∗ soft ( t ) p soft ( t ) , determines the dynamical evolution of the soft sector. One finds that the characteristicmomentum scale p soft ( t ) continues to increase, while the effective temperature T ∗ soft ( t ) ofthe soft sector drops p soft ( t ) ∼ α s f / Q ( Qt ) / , (59) T ∗ soft ( t ) ∼ α − / s f / Q ( Qt ) − / . (60)Eventually, at a time Qt ∼ f − / α − s the characteristic momentum scale p soft ( t ) becomescomparable to the effective temperature T ∗ soft ( t ), indicating that the phase-space densitiesof soft particles f ( p soft ( t )) ∼ e soft ∼ f / e hard of the energy e hard of the hard particles has been transferred to the soft thermal bath via direct radiation. In addition to directly radiating soft gluons with p (cid:46) p soft ( t ),the hard modes can transfer energy to soft sector via multiple successive branchings. Al-though soft branchings with min( z, − z ) (cid:28) z ∼ / p split ( t ) ∼ α s ˆ q ( t ) t , (61)where the probability t dΓ( p split ( t )) / d z to undergo a quasi-democratic splitting with z ∼ / p split ( t ) is the momentum of the most energetic particles that can be
18 S. Schlichting and D. Teaney topped over the total lifetime t of the system (31). Clearly at early times p split ( t ) (cid:28) p soft ( t ),and quasi-democractic branchings contribute very little to the overall energy transfer. How-ever, at a time of order Qt ∼ f − / α − s , p split ( t ) becomes of order of p soft ( t ), and multiplesuccessive branchings begin to dominate the energy transfer to the soft thermal medium.Hence the last stage of the thermalization process is analogous to a highly energetic jetloosing energy to the QGP, highlighting an important connection between jet quenchingand thermalization.In the final stages the soft bath is equilibrated, and ˆ q ( t ) and e soft ( t ) are determined bytheir equilibrium values at temperature T soft ( t ), wich depends on timeˆ q ( t ) ∼ ˆ q soft ( t ) ∼ α s ( T soft ( t )) , e soft ( t ) ∼ ( T soft ( t )) . (62)To determine the rate of energy transfer, we need to compute the energy radiated upto the momentum scale p split ( t ), which will then have time enough to undergo successivebranchings in the bath. Using (53) with the LPM estimate for dΓ / d z from (40), the transferof energy from hard to soft modes is of order e soft ( t ) ∼ (cid:90) t dt (cid:90) p split ( t ) p E p ∂f ( t, | p | ) ∂t ∼ α s e hard (cid:115) ˆ q ( t ) t Q (cid:115) p split ( t ) Q , (63)yielding with (61) the estimate e soft ( t ) ∼ e hard p split ( t ) Q . (64)The transfer of energy ends when the thermal medium has entirely absorbed the energyof the hard partons e soft ( t ) ∼ e hard which occurs when p split ( t ) ∼ Q . Self-consistentlydetermining the time evolution of the scales according to eq. (64) and (62), we find p split ( t ) ∼ α s f Q ( Qt ) and T soft ( t ) ∼ α s f Q ( Qt ) , thus, at a time of order t eq ∼ α − s f − / Q − ∼ α − s T − (cid:114) QT , (65)the temperature of the soft thermal bath T soft ( t ) becomes of the order of the final equilibriumtemperature T ∼ f / Q . In contrast to the overoccupied case, the equilibration time ofan underoccupied system t eq ∼ α − s T − (cid:112) Q/T is parameterically larger than the near-equilibrium relaxation rate ∼ α − s T . Notably, the additional dependence on the ratio ofmomentum scales (cid:112) Q/T implies that excitations with different energies Q equilibrate ondifferent time scales.Beyond the level of parametric estimates (36) a more quantitative description of theinverse energy cascade has been put forward already in the original bottom-up paper (6); theconnections to wave turbulence were established in subsequent works (32, 57) in the contextof jet quenching. Within an inertial range of momenta T soft ( t ) (cid:28) | p | (cid:28) Q the dynamicsis dominated by successive branchings, as described by an effective kinetic equation of theform ∂∂t f ( t, | p | ) (cid:39) (cid:90) dz (cid:20) z − dΓ LPM inel ( p /z )d z f (cid:16) t, p z (cid:17) −
12 dΓ
LPM inel ( p )d z f ( t, p ) (cid:21) , (66)By exploiting the symmetry dΓ LPM inel ( p , z ) = dΓ LPM inel (cid:16) p , − z (cid:17) , and the approximate scaleinvariance of the splitting rates dΓ LPM inel ( p /z, z ) (cid:39) √ z dΓ LPM inel ( p , z ), the collision integral in ••
LPM inel ( p )d z f ( t, p ) (cid:21) , (66)By exploiting the symmetry dΓ LPM inel ( p , z ) = dΓ LPM inel (cid:16) p , − z (cid:17) , and the approximate scaleinvariance of the splitting rates dΓ LPM inel ( p /z, z ) (cid:39) √ z dΓ LPM inel ( p , z ), the collision integral in •• The First fm/c of Heavy-Ion Collisions 19
66) an be transformed into ∂∂t f ( t, p ) (cid:39) (cid:90) dz dΓ LPM inel ( p )d z (cid:104) z − / f (cid:16) t, p z (cid:17) − zf ( t, p ) (cid:105) . (67)Eq. (67) admits stationary solutions of the Kolmogorov-Zakharov form f KZ (cid:16) t, p (cid:17) = f ∗ (cid:18) Q | p | (cid:19) κ , (68)with a universal spectral index κ = 7 / f ∗ . One finds that– in analogy to the Kolmogorov-Zakharov spectra of weak wave turbulence – the solutionis associated with scale independent energy flux, meaning that the energy lost by modesabove a scale Λ ddt e hard ( t ) (cid:39) (cid:90) ∞ Λ πp dp E p ∂∂t f ( t, p ) , (69)is independent of Λ. This property reflects the transport of energy from hard modes (Λ ∼ Q )all the way to the soft thermal bath (Λ ∼ T soft ( t )) via successive branchings. By exploitingthe scale invariance of the collision integral dΓ LPM inel ( p , z ) (cid:39) (cid:112) Q/ | p | dΓ LPM inel ( Q, z ), the energyflux in eq. (69) can be evaluated by using eq. (67) and eq. (68) to evaluate ∂ t f , and bytaking the limit where the spectral exponent approaches the Kolmogorov-Zakharov solutionfrom above (50, 58), κ (cid:38) /
2. This yields ddt e hard ( t ) (cid:39) − (4 π ) Q f ∗ γ g , γ g = Q − (cid:90) dz dΓ LPM inel ( Q, z )d z z log(1 /z ) . (70)While the inverse energy cascade is ultimately responsible for transferring the energy ofhard particles to the soft bath, coincidentally the properties of the QCD splitting functionsare such that a single emission is sufficient to create the turbulent spectrum eq. (68) withinthe inertial range of momenta T soft ( t ) (cid:28) | p | (cid:28) Q (6, 32, 57, 59). Based on this peculiarproperty, it is then also possible to estimate the amount of energy injected into the cascade(corresponding to the non-universal amplitude f ∗ ) and calculate the energy transfer to thethermal bath as discussed in detail in (6, 59). We also note that numerical studies of thethermalization of underoccupied systems performed in (49) confirm the basic picture of thethermalization mechanism illustrated in Fig. 3 and provide additional information on thethermalization time. So far we have discussed the thermalization process for statistically isotropic plasmas. Whenthe distribution is anisotropic, a quantitative analysis of the evolution becomes significantlymore complicated due to the presence of plasma instabilities (60, 61). Once the phase spacedistribution has an order one anisotropy, instabilities qualitatively change the screeningmechanisms in the plasma, and significantly complicate the calculation of radiation ratesand the relaxation to equilibrium (62). How precisely plasma instabilities modify the ther-malization process in over-occupied and under-occupied systems has not been fully clarified,although a number of proposals exist (63, 64). However, it is known that such instabilitiesare much less important than in QED plasmas since the non-linear non-abelian characterof the field equations ultimately limits the growth of the instability (65, 66). While for
20 S. Schlichting and D. Teaney veroccupied systems first-principles studies including the dynamics of instabilities couldbe performed with classical-statistical field simulations, these simulations are technicallychallenging, and most studies in this context have focused on the growth of instabilities atvery early times. Since the situation remains somewhat inconclusive – especially with re-gards to underoccupied systems where classical-statistical simulations are inapplicable – wewill ignore the effects of plasma instabilities throughout the remainder of this section, andonly comment on selected results in our outline of the original bottom-up picture. Currentimplementations in kinetic theory also have ignored plasma instabilities to date (67).As discussed qualitatively in Sect. 2, the first and second/third stages of the bottom-upscenario are characteristic of overoccuppied and underoccupied systems respectively. In allstages the presence of the longitudinal expansion modifies the rates discussed in Sect. 4.1and Sect. 4.2 for static systems, without changing the overall picture. In Fig. 4 we show asimulation result of Kurkela and Zhu of the original bottom-up scenario (67). The simulationuses the ’t Hooft coupling λ = 4 πα s N c (and thus a “realistic” coupling is λ (cid:39)
10 or more ),and starts from a CGC motivated initial condition characterized by1 ν g dNd x ⊥ dy = 0 . Q s λ , (cid:113) (cid:104) p T (cid:105) = 1 . Q s , (71)treating screening with one overall mass m given by eq. (9). The pressure anisotropy isdefined from the stress tensor P T /P L ≡ ( T xx + T yy ) / (2 T zz ), while the occupancy in unitsof λ − is λ (cid:104) pf p (cid:105)(cid:104) p (cid:105) = λ (cid:82) p | p | f p (cid:82) p | p | f p , (72)which in equilibrium reaches 0 . λ , indicated by the crosses in Fig. 4.The numerical simulations confirm the three stage picture of bottom up thermalization:in stage one the the anisotropy grows and the occupancy decreases; in stage two the occu-pancy decreases and the anisotropy is stabilized; and finally in stage three the anisotropyapproaches unity and the energy of the system is thermalized. In the weak coupling limit( λ (cid:39) .
5) the three different stages are clearly visible, whereas for more realistic couplingstrength ( λ (cid:39)
10) the distinctions between the different stages becomes increasingly washedout. We will describe each stage more completely below using the results of Sect. 4.1 andSect. 4.2.To analyze the first over-occupied stage in the expanding case we examine the Boltzmannequation with an elastic scattering (cid:18) ∂∂τ − p z τ ∂∂p z (cid:19) f ( τ, p z , p ⊥ ) = ˆ q ∂ f∂p z . (73)Here the free streaming term on the l.h.s. stems from the expansion of the system, andmakes the momentum distribution increasingly anisotropic (68). On the r.h.s. is the Fokker-Planck operator discussed in Sect. 3, but here we have kept only the most relevant termwhich competes with the expansion and broadens the momentum distribution. Eq. (73)admits a scaling solution of the form f ( τ, p z , p T ) = 1 α ( Q s τ ) / f S (cid:18) p T Q s , p z ( Q s τ ) / Q s (cid:19) , (74) In terms of macroscopic properties, the shear viscosity of the simulation is η/s (cid:39) .
62 for λ = 10. ••
62 for λ = 10. •• The First fm/c of Heavy-Ion Collisions 21 ➀➁➂ λ =0.5 λ =1 λ =10 λ =5 P r e ss u r e an i s o t r op y : P T / P L Phase-space occupancy: λ〈 pf 〉 / 〈 p 〉 Figure 4
Kinetic theory simulation of the non-equilibrium evolution of the pressure anisotropy andphase-space occupancy (see eq. (72) and surrounding text) for a pure Yang-Mills plasma in theoriginal bottom-up scenario (67). Here λ = 4 πα s N c is the coupling, and the black crosses indicateequilibrium value. The three arrows and associated circled numbers indicate the three stages ofbottom-up. provided one uses the by now familiar estimate for ˆ q dominated by the hard modes, ˆ q ∼ α (cid:82) p f p (1 + f p ). This scaling solution, which features a decreasing occupancy and anincreasing anisotropy, is clearly seen in the numerical simulations of (67) at least for thesmallest couplings.The first over-occupied stage of the bottom up scenario has also been addressed indetail within classical-statistical simulations (48, 28, 54). It was found that the phasespace distribution in the classical simulations reaches the universal scaling form of eq. (74),reflecting the NTFP discussed in Sect. 4.1. In these simulations the effects of plasmainstabilities are clearly observed at early times during the approach to the NTFP attractor,but do not appear to significantly affect the longitudinal momentum broadening in thescaling regime, such that (cid:104) p z (cid:105) ∼ Q s ( Q s τ ) − / decreases at late times as in the originalbottom scenario. It remains an open question why plasma instabilities do not seem to playa more important role during the first phase of bottom up.From the scaling solution in eq. (74), we see that the first phase ends at a time Q s τ ∼ α − / s , and after this point the system in is an under-occupied non-equilibrium state. Theestimates and physics for the thermalization of such states described in Sect. 4.2 can beadapted to the expanding case by recognizing that hard modes are essentially free streaming,and thus the energy and number densities of these modes are continually decreasing, so thatthe energy and number per rapidity ( τ e and τ n respectively) remains fixed: τ e hard ( τ ) = Q s α s , (75) τ n hard ( τ ) = Q s α s . (76)
22 S. Schlichting and D. Teaney sing the estimate ˆ q ( τ ) ∼ α s (cid:90) p f p (1 + f p ) ∼ α s n hard ( τ ) , (77)one finds that because of the expansion the soft scale p soft ( τ ) remains constant in time p ( τ ) ∼ ˆ q ( τ ) τ ∼ α s Q s , (78)as opposed to increasing as it does in the non-expanding case. Thus, the pressure anisotropyin the second phase is constant and large as seen in Fig. 4. Eq. (58) for the energy densityproduced by direction radiation by the bath into the soft modes remains valid e soft ( τ ) ∼ α s τ e hard ( τ ) (cid:115) ˆ q ( τ ) Q s (cid:115) p soft ( τ ) Q s , (58)but now e hard ( τ ) and ˆ q ( τ ) are functions of time. Qualitatively, Eq. 58 will hold even ifplasma instabilities are present, but ˆ q will deviate from the estimate in eq. (77), which isbased upon elastic scattering by the hard modes. However, because the plasma instabilitiesare bounded they will not radically change the picture. The second phase of bottom-up endswhen e soft ( τ ) ∼ p ( τ ) and the soft bath has thermalized. Equating these two expressionone finds that the second phase ends at a time of order Q s τ ∼ α − / s .Finally, we analyze the last phase of bottom-up. Here again the physics is identical tothe inverse energy cascade discussed in detail in Sect. 4.2.2 for the static system. Eq. (61)for the splitting (or stopping) momentum p split ( τ ) = α s ˆ q ( τ ) τ , and eq. (64) for e soft ( τ ) areunchanged e soft ( τ ) ∼ e hard ( τ ) p split ( τ ) Q s , (64)provided the free streaming result for e hard ( τ ) is used. Again, plasma instabilities may mod-ify our estimate for ˆ q ( τ ), but this will not change the overall picture. The system is com-pletely thermalized when e soft ( τ ) becomes comparable to e hard ( τ ), τ e soft ( τ ) ∼ τ e hard ( τ ) ∼ Q s /α . Using the fact that ˆ q is determined by e soft in equilibrium, ˆ q ( τ ) ∼ α s e / ( τ ), onereadily establishes that the system thermalizes at Q s τ ∼ α − / s . (79)We hope that it is evident that the overall picture of bottom-up is quite robust. Ulti-mately this picture follows from a hard scale Q s , kinematics, and generic features of collinearradiation. These features tend to fill up a soft sector first, which then causes a cascade ofthe energy of the system to the IR. Indeed, an extensive analysis of thermalization whenplasma instabilities are present finds many of the same qualitative features of bottom-upwith somewhat modified exponents (63).
5. Simulations of early time dynamics and heavy-ion phenomenology5.1. Approach to hydrodynamics
We now turn to simulations of the early time dynamics and the approach to equilibrium inhigh-energy heavy ion collisions. Here we will focus on the eventual approach towards localthermal equilibrium, and determine when the evolution can be described with relativisticviscous fluid dynamics. ••
We now turn to simulations of the early time dynamics and the approach to equilibrium inhigh-energy heavy ion collisions. Here we will focus on the eventual approach towards localthermal equilibrium, and determine when the evolution can be described with relativisticviscous fluid dynamics. •• The First fm/c of Heavy-Ion Collisions 23 iscous fluid dynamics describes the macroscopic evolution of the energy-momentumtensor T µν , based on an expansion around local thermal equilibrium which is controlledby the Knudsen number Kn θ ∼ τ micro / T macro , and the proximity to the equilibriumstate, which can be quantified by the non-equilibrium corrections to the stress tensor, T µν non − eq /T µν eq . At early times τ ∼ /Q s the longitudinal pressure is much smaller thanthe transverse pressure P L (cid:28) P T and hydrodynamics does not apply. Consequently, thekey question is to understand how T µν then evolves towards local thermal equilibrium wherethe longitudinal and transverse pressures are equal P L = P T .Neglecting potential problems related to plasma instabilities, the non-equilibrium evolu-tion of macroscopic quantities such as T µν can be calculated based on numerical simulationsof the effective kinetic theory. Numerical simulation based on QCD kinetic theory werepioneered in (69, 70); the first complete leading order study for a homogeneous purely glu-onic plasma was performed in (67) and subsequently extended to inhomogeneous plasmas(71, 30, 72) as well as homogeneous plasmas of quarks and gluons (73, 74). Kinetic theorysimulations shown in Fig. 5(a) indicate that for realistic coupling strength α s (cid:38) .
1, theevolution of the energy-momentum tensor towards equilibrium is to a good approximationcontrolled by a single time scale τ eq R , corresponding to the equilibrium relaxation rate τ eq R ( τ ) = 4 πη/sT Id ( τ ) , (80)where η/s ∝ λ is the shear-viscosity to entropy density ratio, and T Id ( τ ) = ∝ τ − / de-notes the temperature of the late-time equilibrium system. Even though an extrapolationto sizeable coupling strength is required to make contact with heavy-ion phenomenology,the dependence on α s is surprisingly weak once τ is measured in units of τ eq R . When com-paring the results for the non-equilibrium evolution of the energy-momentum tensor T µν in kinetic theory with the asymptotic behavior in viscous hydrodynamics, one concludesthat a fluid dynamic description becomes applicable on time scales τ hydro ≈ τ eq R ( τ ). Forphenomenological purposes the coupling constant λ can be traded for η/s ∝ λ yielding thefollowing estimate τ hydro ≈ . (cid:18) π ( η/s )2 (cid:19) / (cid:18) (cid:104) τ s (cid:105) . (cid:19) − / (cid:16) ν eff (cid:17) / , (81)where (cid:104) τ s (cid:105) denotes the entropy density per unity rapidity. τ s is directly related tothe charged particle multiplicity dN ch /dη , and thereby constrained to be approximately ≈ . for central Pb+Pb collisions at LHC energies (71). Since the discussion so farignores the effects of spatial gradients, both in transverse space and longitudinal rapidity,the estimate in (81) should be understood as a lower bound.Interestingly one finds that viscous hydrodynamics starts to describe the evolution ofthe energy-momentum tensor in a regime where both the Knudsen number Kn θ ≈ τ eq R /τ andthe proximity to equilibrium as measured by 1 − P L /P T are of order unity, indicating thatthe system is still significantly out-of-equilibrium. Even though this behavior appears to bequite surprising, it is by no means unique to a weakly coupled non-equilibrium description,and similar observations have been reported much earlier in the context of strongly coupled T macro is a typical macro timescale, which can be estimated from the inverse of expansion scalar( ∇ · u ) ≡ T − of the fluid. For a Bjorken expansion T macro = τ .
24 S. Schlichting and D. Teaney . . . .
81 0 1 2 3 4 5 ¯ T µ ν / ¯ T µ ν i d . τT id . / (4 πη/s ) λ = 10 ( η/s ≈ . ) λ = 15 ( η/s ≈ . ) λ = 20 ( η/s ≈ . ) λ = 25 ( η/s ≈ . ) e/ P T P L . . . . τ hydro τ chem τ therm η/s = 0 . e [ G e V / f m ] τ [fm] totalgluonsfermions Figure 5 (left) Non-equilibrium evolution of the different components of the average energy-momentumtensor T µν = diag( e, P T , P T , P L ) compared to viscous fluid dynamics (30, 72). (right) Evolutionof the overall energy density e and the energy density carried by quarks and gluons e g/q (73). gauge theories (75). It has become common to distinguish the time when hydrodynamicsbecomes applicable τ hydro (the so called “hydrodynamization” time) from the time τ eq when the pressure anisotropy is small. Due to the rapid longitudinal expansion, the actualapproach towards local pressure isotropy occurs only on much larger time scales τ eq (cid:29) τ hydro . Hence the great success of hydrodynamic descriptions of the QGP does not appearto derive from the fact that the system is particularly close to equilibrium throughoutmost of its space-time evolution, but is rather due to fact that the range of applicabilityof viscous relativistic fluid appears to be larger than originally anticipated. Notably, theseobservations have triggered a large number of theoretical studies to further investigate andpossibly extend the range of applicability of viscous fluid dynamics (76, 77, 78). However,a detailed discussion of these topics is beyond the scope of this review. So far most theoretical studies of the early non-equilibrium dynamics have focused onthe kinetic equilibration of gluons, while neglecting dynamical fermions in the analysis.However, on a conceptual level it is equally important to understand the transition from aninitial state, which is believed to be highly gluon dominated, towards chemical equilibriumwhere a significant fraction of the energy density is carried by quark degrees of freedom.We note from a phenomenological point of view the chemical composition of the plasmaat early times, may have also have interesting consequences, e.g. relating to the questionsconcerning the chemical equilibration of strange quarks and heavy flavors or the electro-magnetic response of the QGP at very early times after the collision. Even though acomplete picture of chemical equilibration along the lines of our discussion in Sec. 4.2 is yetto be established, interesting first results have been reported in the literature. We brieflydiscuss these results below.Classical-statistical simulations of quark production at very early times have been pi-oneered in (79) demonstrating that at realistic coupling strength a significant number ofquark anti-quark pairs can be produced in the initial (semi-) hard scattering and in the ••
81 0 1 2 3 4 5 ¯ T µ ν / ¯ T µ ν i d . τT id . / (4 πη/s ) λ = 10 ( η/s ≈ . ) λ = 15 ( η/s ≈ . ) λ = 20 ( η/s ≈ . ) λ = 25 ( η/s ≈ . ) e/ P T P L . . . . τ hydro τ chem τ therm η/s = 0 . e [ G e V / f m ] τ [fm] totalgluonsfermions Figure 5 (left) Non-equilibrium evolution of the different components of the average energy-momentumtensor T µν = diag( e, P T , P T , P L ) compared to viscous fluid dynamics (30, 72). (right) Evolutionof the overall energy density e and the energy density carried by quarks and gluons e g/q (73). gauge theories (75). It has become common to distinguish the time when hydrodynamicsbecomes applicable τ hydro (the so called “hydrodynamization” time) from the time τ eq when the pressure anisotropy is small. Due to the rapid longitudinal expansion, the actualapproach towards local pressure isotropy occurs only on much larger time scales τ eq (cid:29) τ hydro . Hence the great success of hydrodynamic descriptions of the QGP does not appearto derive from the fact that the system is particularly close to equilibrium throughoutmost of its space-time evolution, but is rather due to fact that the range of applicabilityof viscous relativistic fluid appears to be larger than originally anticipated. Notably, theseobservations have triggered a large number of theoretical studies to further investigate andpossibly extend the range of applicability of viscous fluid dynamics (76, 77, 78). However,a detailed discussion of these topics is beyond the scope of this review. So far most theoretical studies of the early non-equilibrium dynamics have focused onthe kinetic equilibration of gluons, while neglecting dynamical fermions in the analysis.However, on a conceptual level it is equally important to understand the transition from aninitial state, which is believed to be highly gluon dominated, towards chemical equilibriumwhere a significant fraction of the energy density is carried by quark degrees of freedom.We note from a phenomenological point of view the chemical composition of the plasmaat early times, may have also have interesting consequences, e.g. relating to the questionsconcerning the chemical equilibration of strange quarks and heavy flavors or the electro-magnetic response of the QGP at very early times after the collision. Even though acomplete picture of chemical equilibration along the lines of our discussion in Sec. 4.2 is yetto be established, interesting first results have been reported in the literature. We brieflydiscuss these results below.Classical-statistical simulations of quark production at very early times have been pi-oneered in (79) demonstrating that at realistic coupling strength a significant number ofquark anti-quark pairs can be produced in the initial (semi-) hard scattering and in the •• The First fm/c of Heavy-Ion Collisions 25 resence of the strong color fields at very early times. Subsequent studies have refinedthe lattice approach (80, 81) and further elaborated on quark production in over-occupiedsystems (82). However, as classical-statistical simulations involving dynamical fermions aresignificantly more complicated, studies are yet to reach the same level of sophistication ofanalogous pure gauge theory simulations.Quark production during the final radiative break-up stage of the bottom up scenario,has been investigated in the context of jet quenching (59), where it was pointed out thatthe turbulent nature of the inverse energy cascade ultimately determines the quark/gluonratio from a local balance of the g → q ¯ q and q → qg processes. However, within the inertialrange of the cascade T soft (cid:28) p (cid:28) p split the fraction of energy carried by quarks and anti-quarks e q /e g (cid:39) . × N f (for three colors) is small compared to the equilibrium ratio e q /e g (cid:39) . × N f , indicating that elastic processes, which are operative at the scales ofthe soft thermal medium also play a pivotal role in the chemical equilibration process.The first numerical study implementing all relevant leading order processes of bottomup was performed in (73, 74), indicating that as shown in Fig. 5 the approach to viscousfluid dynamic behavior (discussed in Sect. 5.1) occurs before chemical equilibration of theQGP. A complete leading order analysis of the chemical equilibration mechanism (along thelines of Sect. 4.2) has not yet been given, and should explain these first numerical resultsand provide guidance to phenomenology.Notably, the inclusion of dynamical quarks also represent an important step towardscalculations of pre-equilibrium photon and dilepton production, and in addressing questionsrelated to the chemical/kinetic equilibration of heavy flavors. While first progress in thisdirection has been reported in (83) by analyzing a subset of leading order processes, acomplete leading order study has not been performed to date. So far we have discussed the microscopic dynamics of the local equilibration process, ne-glecting the effects of spatial gradients on small scales ∼ cτ hydro . However, as discussed inSect. 2 the inclusion of small scale fluctuations ∼ R p is a necessary ingredient for a realisticevent-by-event description, since such gradients will lead to the development of “pre-flow”,a pre-cursor to the late stage hydrodynamic flow which starts to build up already during thepre-equilibrium phase. The kinetic theory should evolve these fluctuations and smoothlyasymptote to hydrodynamics at late times τ ∼ τ hydro .A recent extension of the bottom up scenario accounts for small scale fluctuations byexplicitly including spatially inhomogeneous fluctuations of the phase space density intothe kinetic description (30, 72, 71). By choosing a representative form for the phase-spacedistribution to model the initial fluctuations of the stress tensor δT µν ( τ , x ) around a localaverage ¯ T µν x ( τ ) at a point x , the pre-equilibrium evolution of the energy-momentum tensorcan then be calculated as T µν ( τ, x ) (cid:39) ¯ T µν x ( τ ) (cid:124) (cid:123)(cid:122) (cid:125) non-eq. evolution of(local) avg. background + (cid:90) (cid:12) d x G µναβ ( τ, τ , x , x ) δT αβ ( τ , x ) (cid:124) (cid:123)(cid:122) (cid:125) non-eq. evolution of local fluctuations of the stress tensor , (82)which is shown schematically in Fig. 6. Here ¯ T µν x ( τ ) describes the pre-equilibrium evolutionof the average energy-momentum tensor and is described by Fig. 5, while G µναβ describesthe linear response to initial fluctuations in the thermalizing plasma (30). Since causality
26 S. Schlichting and D. Teaney τ = τ Hydro τ = τ EKT ~ s Initial stateHydro K i ne t i c Theo r y Figure 6
Schematic of the transverse energy density profiles at very early times ( τ = τ EKT ≈ . /c ) andafter the first fm /c of pre-equilibrium evolution done with kinetic theory ( τ = τ hydro ). At a time τ hydro the constitutive equations are approximately satisfied (see Fig. 7). restricts the contributions to the fluctuations at x , one only needs to integrate the responseover the causal circle (cid:12) indicated by the circle at τ = τ EKT in Fig. 6. The relevant responsefunctions ¯ T µν ( τ ) and G µναβ , can be calculated once and for all in kinetic theory, and packagedinto a useful “pre-flow” computer code which encapsulates the thermalization process (30).The linear response formalism of eq. (82) can be seen as a systematic extension of ear-lier studies (84), recognizing universal patterns in the pre-equilibrium evolution of the longwave-length components of the energy-momentum tensor. Short wave-length fluctuations (cid:28) cτ Hydro are efficiently damped during the pre-equilibrium phase, leading to an effectivecoarse graining of the spatial profile of the energy-momentum tensor shown schematicallyin Fig. 6. Then viscous corrections to the energy-momentum tensor are reasonably wellapproximated by the Navier-Stokes constitutive relations at the time τ init . when hydrody-namics is initialized. This is shown in Fig. 7 (left), which uses eq. (82) to determine thestress at a time τ init . . Long wave-length fluctuations of the initial energy density deter-mine the pre-flow which develops during thermalization process, and can be reasonablyapproximated as T τi ( τ, x ) ≈ − ( τ − τ )2 (cid:18) ¯ T ττ x ( τ )¯ T ττ x ( τ ) (cid:19) ∂ i ¯ T ττ ( τ , x ) . (83)Nevertheless, the results of (30, 72) also demonstrate that a genuine non-equilibrium de-scription is necessary account for the entropy production during the pre-equilibrium phase.Since the subsequent hydrodynamic expansion approximately conserves the overall entropy,this factor two to three increase in entropy during the pre-equilibrium phase is importantin relating properties of the initial state to experimentally observed charged particle multi-plicities. ••
Schematic of the transverse energy density profiles at very early times ( τ = τ EKT ≈ . /c ) andafter the first fm /c of pre-equilibrium evolution done with kinetic theory ( τ = τ hydro ). At a time τ hydro the constitutive equations are approximately satisfied (see Fig. 7). restricts the contributions to the fluctuations at x , one only needs to integrate the responseover the causal circle (cid:12) indicated by the circle at τ = τ EKT in Fig. 6. The relevant responsefunctions ¯ T µν ( τ ) and G µναβ , can be calculated once and for all in kinetic theory, and packagedinto a useful “pre-flow” computer code which encapsulates the thermalization process (30).The linear response formalism of eq. (82) can be seen as a systematic extension of ear-lier studies (84), recognizing universal patterns in the pre-equilibrium evolution of the longwave-length components of the energy-momentum tensor. Short wave-length fluctuations (cid:28) cτ Hydro are efficiently damped during the pre-equilibrium phase, leading to an effectivecoarse graining of the spatial profile of the energy-momentum tensor shown schematicallyin Fig. 6. Then viscous corrections to the energy-momentum tensor are reasonably wellapproximated by the Navier-Stokes constitutive relations at the time τ init . when hydrody-namics is initialized. This is shown in Fig. 7 (left), which uses eq. (82) to determine thestress at a time τ init . . Long wave-length fluctuations of the initial energy density deter-mine the pre-flow which develops during thermalization process, and can be reasonablyapproximated as T τi ( τ, x ) ≈ − ( τ − τ )2 (cid:18) ¯ T ττ x ( τ )¯ T ττ x ( τ ) (cid:19) ∂ i ¯ T ττ ( τ , x ) . (83)Nevertheless, the results of (30, 72) also demonstrate that a genuine non-equilibrium de-scription is necessary account for the entropy production during the pre-equilibrium phase.Since the subsequent hydrodynamic expansion approximately conserves the overall entropy,this factor two to three increase in entropy during the pre-equilibrium phase is importantin relating properties of the initial state to experimentally observed charged particle multi-plicities. •• The First fm/c of Heavy-Ion Collisions 27 − − ( π xx + π yy ) G e V / f m x fm τ init . = 0 . τ init . = 1 . τ init . = 1 . − − . . .
01 0 . h P T ih e ih P L ih e i / τ (fm) τ EKT τ hydro Figure 7 (left) Spatial profiles of the non-equilibrium shear-stress tensor (Π xx + Π yy ) (the solid lines)compared to the Navier-Stokes hydrodynamics limit (the dashed lines) after an evolution of τ init . in the kinetic theory (see eq. (82)). (right) Proof of principle calculation combining differenttheoretical descriptions to calculate the evolution of energy density e , P T , P L in a single Pb+Pbevent (30, 72). Finally, by combining the classical-statistical field simulations at early times, the kineticsimulations at intermediate times, and the hydrodynamics simulations at late times, aconsistent space-time description of the energy-momentum tensor can be obtained on anevent-by-event basis. This is illustrated in Fig. 7 (right) which shows the evolution oflongitudinal and transverse pressures in a single Pb+Pb event. In this simulation the firststage up to τ EKT is treated in the classical IP-Glasma model (see Sect. 2); the second stageup to τ hydro is treated with QCD kinetics following the outlines of bottom-up; and thefinal phase is treated with hydrodynamics. The different theoretical descriptions overlapproviding a complete picture of the event.
6. Outlook and small systems
We have reviewed the weak coupling description of the thermalization process of the QGPduring the first fm/c of high-energy heavy-ion collisions, by dividing the out-of-equilibriumdynamics of non-abelian gauge theories into two broad classes – an over-occupied limitdiscussed in Sect. 4.1, and an under-occupied limit discussed in Sect. 4.2. Strikingly, thethermalization process in each of these limits exhibits generic scaling features which onewould like to observe experimentally.Indeed, much of the current interest in the equilibration process is driven by excitingnew data on the small systems created in proton-proton ( p + p ) and proton-nucleus ( p + A )collisions, which show evidence for a transition towards a hydrodynamic regime in nucleus-nucleus A + A collisions. A more complete review of the experimental data is given in theliterature (85, 2). In the larger A + A system, the approach to hydrodynamics has beenlargely understood and quantified within the bottom-up scenario (see for example eq. (81) ofSect. 5), and the physics of the pre-equilibrium stage has been packaged into a useful “pre-flow” computer code that can be used to simulate heavy ion events (see Sect. 5.3). However,as the system size gets smaller, additional scales, such as the transverse radius R , play an
28 S. Schlichting and D. Teaney ncreasingly important role and truncate the thermalization process. Nevertheless, one canuse the bottom-up framework to estimate when hydrodynamics becomes applicable as afunction of the multiplicity produced in the collision (30). Substituting τ s = dS/dy/πR in eq. (81) we find τ hydro R (cid:39) (cid:18) dN ch /dy (cid:19) − / (cid:18) πη/s (cid:19) / (cid:18) S/N ch (cid:19) − / (cid:16) ν eff (cid:17) / . (84)Since we expect the bottom up analysis will be strongly modified for τ hydro /R ∼
1, a chargedparticle multiplicity of order dN ch /dy ∼
70 should demarcate the transition to a fullyequilibrated regime. So far a detailed understanding, both parametrically and numericallyof the transition regime has not been given, though important first steps have been taken (29,86, 87). In small systems there are by now many experimental tools, (such as e.g. the hadronchemistry (73) or the system size dependence of the harmonic flow (88, 89)) which can beused to clarify the kinetics of high energy QCD and to guide theory. Further as emphasizedin Sect. 4.2.2, studies of the energy loss of jets, both in small and large systems, can informthe study of thermalization of QCD plasmas. We therefore anticipate that, through acombination of phenomenology, formal theory, experiment, and simulation, the communitywill analyze the transition from cold QCD to the hot QGP in detail, and, more generally,clarify the out-of-equilibrium behavior of non-abelian gauge theories.
DISCLOSURE STATEMENT
The authors are not aware of any affiliations, memberships, funding, or financial holdingsthat might be perceived as affecting the objectivity of this review.
ACKNOWLEDGMENTS
We gratefully acknowledge helpful discussions with Peter Arnold, Juergen Berges, AleksiKurkela, Aleksas Mazeliauskas, Jean-Francois Paquet, and Raju Venugopalan. This workwas supported in part by the U.S. Department of Energy, Office of Science, Office of NuclearPhysics under Award Numbers DE-FG02-88ER40388 and by the Deutsche Forschungsge-meinschaft (DFG, German Research Foundation) Project number 315477589 TRR 211.
LITERATURE CITED
1. Heinz U, Snellings R.
Ann. Rev. Nucl. Part. Sci.
Int. J. Mod. Phys.
E25:1630002 (2016)3. Nagle JL, Zajc WA.
Ann. Rev. Nucl. Part. Sci.
Acta Phys. Polon.
B47:2581 (2016)5. Keegan L, et al.
JHEP
Phys. Lett.
B502:51 (2001)7. Bjorken JD.
Phys. Rev.
D27:140 (1983)8. Gelis F, Schenke B.
Ann. Rev. Nucl. Part. Sci.
Phys. Rev.
D49:2233 (1994)10. Kovner A, McLerran LD, Weigert H.
Phys. Rev.
D52:3809 (1995)11. Kovner A, McLerran LD, Weigert H.
Phys. Rev.
D52:6231 (1995)12. Lappi T, McLerran L.
Nucl. Phys.
A772:200 (2006)13. Mazeliauskas A, Teaney D.
Phys. Rev.
C91:044902 (2015) ••
C91:044902 (2015) •• The First fm/c of Heavy-Ion Collisions 29
4. Iancu E, Venugopalan R. 2003. In
In *Hwa, R.C. (ed.) et al.: Quark gluon plasma* 249-3363
15. Gelis F, Iancu E, Jalilian-Marian J, Venugopalan R.
Ann. Rev. Nucl. Part. Sci.
Phys. Rev. Lett.
Phys. Rev. Lett.
Phys. Rev.
D88:085015 (2013)19. Romatschke P, Venugopalan R.
Phys. Rev. Lett.
Phys. Rev.
D74:045011 (2006)21. Berges J, Schlichting S.
Phys. Rev.
D87:014026 (2013)22. Schenke B, Tribedy P, Venugopalan R.
Phys. Rev. Lett.
Phys. Rev.
C86:034908 (2012)24. Mueller AH, Son DT.
Phys. Lett.
B582:279 (2004)25. Aarts G, Smit J.
Nucl. Phys.
B511:451 (1998)26. Jeon S.
Phys. Rev.
C72:014907 (2005)27. Berges J, Schlichting S, Sexty D.
Phys. Rev.
D86:074006 (2012)28. Berges J, Boguslavski K, Schlichting S, Venugopalan R.
Phys. Rev.
D89:074011 (2014)29. Greif M, et al.
Phys. Rev.
D96:091504 (2017)30. Kurkela A, et al. arXiv:1805.00961 [hep-ph] (2018)31. Arnold PB, Cantrell S, Xiao W.
Phys. Rev.
D81:045017 (2010)32. Blaizot JP, Iancu E, Mehtar-Tani Y.
Phys. Rev. Lett.
JHEP
JHEP
JHEP
JHEP
Phys. Rept.
Course of theoretical physics . London, Pergamon Press;Reading, Mass., Addison-Wesley Pub. Co. (1958)39. Ghiglieri J, Teaney D.
Int. J. Mod. Phys.
E24:1530013 (2015)40. Ghiglieri J, Moore GD, Teaney D.
JHEP
Phys. Rev.
D79:065025 (2009)42. Arnold PB, Xiao W.
Phys. Rev.
D78:125008 (2008)43. Baier R, et al.
Nucl. Phys.
B483:291 (1997)44. Zakharov BG.
JETP Lett.
Phys. Rev.
D25:746 (1982)46. Arnold PB, Dogan C.
Phys. Rev.
D78:065008 (2008)47. Blaizot JP, et al.
Nucl. Phys.
A873:68 (2012)48. Berges J, Boguslavski K, Schlichting S, Venugopalan R.
Phys. Rev.
D89:114007 (2014)49. Kurkela A, Lu E.
Phys. Rev. Lett.
Phys. Rev.
D86:056008 (2012)52. Schlichting S.
Phys. Rev.
D86:065008 (2012)53. Abraao York MC, Kurkela A, Lu E, Moore GD.
Phys. Rev.
D89:074036 (2014)54. Berges J, Mace M, Schlichting S.
Phys. Rev. Lett.
Phys. Rev.
D93:074036 (2016)56. Boguslavski K, Kurkela A, Lappi T, Peuron J.
Phys. Rev.
D98:014006 (2018)57. Blaizot JP, Mehtar-Tani Y.
Annals Phys.
JHEP
Phys. Lett.
B314:118 (1993)61. Romatschke P, Strickland M.
Phys. Rev.
D68:036004 (2003)62. Arnold PB, Lenaghan J, Moore GD.
JHEP
30 S. Schlichting and D. Teaney
3. Kurkela A, Moore GD.
JHEP
JHEP
Phys. Rev. Lett.
Phys. Rev.
D72:054003 (2005)67. Kurkela A, Zhu Y.
Phys. Rev. Lett.
Phys. Lett.
Phys. Rev.
C71:064901 (2005)70. El A, Xu Z, Greiner C.
Nucl. Phys.
A806:287 (2008)71. Keegan L, Kurkela A, Mazeliauskas A, Teaney D.
JHEP
Phys. Rev. Lett.
Rept. Prog. Phys.
Phys. Rev. Lett.
Phys. Rev.
D97:036020 (2018)79. Gelis F, Kajantie K, Lappi T.
Phys. Rev. Lett.
JHEP
Phys. Rev. Lett.
Phys. Rev.
D97:034013 (2018)83. Berges J, Reygers K, Tanji N, Venugopalan R.
Phys. Rev.
C95:054904 (2017)84. Vredevoogd J, Pratt S.
Phys. Rev.
C79:044915 (2009)85. Loizides C.
Nucl. Phys.
A956:200 (2016)86. Borghini N, Feld S, Kersting N.
Eur. Phys. J.
C78:832 (2018)87. Kurkela A, Wiedemann UA, Wu B.
Phys. Lett.
B783:274 (2018)88. Yan L, Ollitrault JY.
Phys. Rev. Lett.
Phys. Lett.
B788:161 (2019) ••