Broad excitations in a 2+1D overoccupied gluon plasma
LLU-TP 21-01
Broad excitations in a 2+1D overoccupied gluon plasma
K. Boguslavski, A. Kurkela, T. Lappi,
3, 4 and J. Peuron
5, 6 Institute for Theoretical Physics, Technische Universit¨at Wien, 1040 Vienna, Austria Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway Department of Physics, University of Jyv¨askyl¨a,P.O. Box 35, 40014 University of Jyv¨askyl¨a, Finland Helsinki Institute of Physics, P.O. Box 64, 00014 University of Helsinki, Finland Dept. of Astronomy and Theoretical Physics, S¨olvegatan 14A, S-223 62 Lund, Sweden European Centre for Theoretical Studies in Nuclear Physics and Related Areas (ECT*)and Fondazione Bruno Kessler, Strada delle Tabarelle 286, I-38123 Villazzano (TN), Italy
Motivated by the initial stages of high-energy heavy-ion collisions, we study excitations of far-from-equilibrium 2+1 dimensional gauge theories using classical-statistical lattice simulations. Weevolve field perturbations over a strongly overoccupied background undergoing self-similar evolution.While in 3+1D the excitations are described by hard-thermal loop theory, their structure in 2+1D isnontrivial and nonperturbative. These nonperturbative interactions lead to broad excitation peaksin spectral and statistical correlation functions. Their width is comparable to the frequency ofsoft excitations, demonstrating the absence of soft quasiparticles in these theories. Our results alsosuggest that excitations at higher momenta are sufficiently long-lived, such that an effective kinetictheory description for 2+1 dimensional Glasma-like systems may exist, but its collision kernel mustbe nonperturbatively determined.
I. INTRODUCTION
In the weak coupling framework, the initial stage of arelativistic heavy-ion collision in the high-energy limit isdominated by nonperturbatively strong boost invariantgluon fields [1]. At the characteristic momentum scale Q s , these strong fields have nonperturbatively large oc-cupation numbers f ∼ /g and are thus effectively clas-sical [2, 3]. The success of a hydrodynamic description ofthe later stages of the evolution indicates that the mat-ter evolves towards sufficient isotropy quickly. One ofthe major open questions in the field is to understandin a controlled theoretical framework the evolution froma classical field-dominated system with approximatelyboost-invariant fields at τ (cid:46) /Q s to a state close tolocal thermal equilibrium [3].Several studies in the past have addressed the evolutionof the initial boost-invariant field during the timescale τ ∼ /Q s , which can be understood in terms of 2+1-dimensional lattice simulations [4–7] or analytical calcu-lations [2, 8–10]. These approaches have been combinedwith hydrodynamic evolution in the later stages, pro-ducing a phenomenologically successful descriptions ofheavy-ion collisions [11–13]. However, these calculationshave not been based on a detailed microscopical descrip-tion of the evolution from boost invariant classical fieldsto hydrodynamics.In practice, the initial boost invariance is broken by dif-ferent mechanisms. At finite collision energy it is brokenby the finite thickness of the nucleus [14–16]. Anothersource of violation are rapidity dependent unstable quan-tum fluctuations [17–21]. After the boost invariance hasbeen broken and the system has become more dilute withan occupation number 1 (cid:28) f (cid:28) /g , the classical fieldscan also be described using kinetic theory [22–25]. Thisstage has been studied extensively using real time lattice simulations [25–31]. More recent studies, which matchkinetic theory to hydrodynamics [32–34], seem to favor abottom-up type isotropization scenario [35]. However, itis unclear if the kinetic theory approach can be used allthe way back to the purely 2+1-dimensional (2+1D) ini-tial stage. This is one of the questions we aim to addressin this paper, which is a continuation of our earlier studyof equal-time correlation functions in the same systemsin Ref. [36].Formulating a kinetic theory in 2+1D faces certainproblematic issues [36]. The kinetic theory relies on aseparation of modes, into soft modes described by classi-cal field theory and hard momentum modes described asparticles that undergo scattering. In order to find lead-ing order expressions of the medium-modified scatteringcross-sections of the hard particles, the dynamics of thesoft sector need to be solved. In three dimensions thedominant medium modification to the soft fields arisesfrom the interaction with the hard modes, and this inter-action can be analytically solved giving rise to the Hard-Loop (HTL) theory [37, 38] and the Debye mass scale m D . However, in two spatial dimensions, the lower di-mensionality puts more weight on infrared modes in mo-mentum integrals, and consequently, the soft and hardscales contribute equally to the Debye screening mass.As a consequence of their self-interactions, already modesat the soft scale m D become nonperturbative, bearingresemblance to the magnetic scale in three dimensions.This means that it is not obvious whether the soft modescan be described using HTL theory. This makes it dif-ficult to find leading order accurate matrix elements forkinetic theory in 2+1D analytically.Nevertheless, our recent studies in 2+1D classical lat-tice gauge theory [36] indicated that not only isotropic3+1D but also 2+1D systems exhibit a self-similar at-tractor behavior, for which we extracted the scaling ex- a r X i v : . [ h e p - ph ] J a n ponents. We also found that the scaling exponents can beunderstood using a simple kinetic theory analysis. Thissuggests that 2+1D systems involve quasiparticle excita-tions, at least at high momenta. Inspired by these results,our aim in this paper is to study the excitation spectrumof the pure 2+1D gauge theory. We will study theorieswith and without a scalar field in the adjoint represen-tation. Especially the former case closely resembles theboost invariant initial state of heavy-ion collisions at highenergy, as the scalar field arises by dimensional reductionfrom 3+1 dimensional pure gauge theory.In practice, the excitation spectrum is studied usingthe linear response framework built in [39] and first ap-plied in [40] to an isotropic 3+1D gluonic system. Asimilar classical-statistical linear response framework hasbeen recently used in scalar theories at self-similar attrac-tors [41, 42] and extends classical-statistical simulationsin thermal equilibrium that use a fluctuation-dissipationrelation explicitly [43, 44]. Our work is a natural exten-sion of our previous study of the excitation spectrum ofisotropic 3+1D gluodynamics at a classical self-similarattractor [40]. There we have observed narrow quasi-particle excitation peaks in the spectrum for all mo-menta and a generalized fluctuation-dissipation relationbetween spectral and statistical correlation functions.While we also find a generalized fluctuation-dissipationrelation, our main result is that the gluonic correlationfunctions in 2+1D systems exhibit a broad excitationpeak of non-Lorentzian shape, in contrast to 3+1D. Infact, we show that the peak width is of the order ofand scales with the mass. We separately confirm in Ap-pendix B that this happens also in classical thermal equi-librium. This indicates that gluonic excitations for mo-menta below the mass are too short-lived to form quasi-particles. Moreover, expressions obtained from HTL the-ory mostly do not provide a good description of our sim-ulation results. This is in line with the breakdown of theHTL theory at the soft scale. The interaction betweensoft excitations is genuinely nonperturbative in 2+1D,and the dynamical description has to take this into ac-count. Our results also suggest that gluonic and scalarexcitations at higher momenta are sufficiently long-lived,such that we may have an effective kinetic descriptionalso in the strongly anisotropic systems which arise in theearly stages of ultrarelativistic heavy ion collisions, albeitwith nonperturbatively detemined collision kernels.We start in Sec. II by introducing the consideredtheories and relevant correlation functions, after whichwe describe our numerical method and introduce themain HTL expressions. Our numerical results for non-equilibrium 2+1D systems are shown in Sec. III. We con-clude in Sec. IV. The analytical forms for the HTL func-tions in 2+1D are derived in Appendix A, and a numeri-cal study in classical thermal equilibrium is discussed inAppendix B. II. THEORETICAL BACKGROUNDA. 3+1D, 2+1D and Glasma-like 2+1D theories
We consider non-Abelian SU( N c ) gauge theories with N c = 2 in d spatial dimensions. Their classical actionreads S Y M [ A ] = − (cid:90) d d +1 x F µνa F aµν , (1)with field strength tensor F aµν = ∂ µ A aν − ∂ ν A aµ + g f abc A bµ A cν , where repeated color indices a = 1 , ..., N c − µ, ν = 0 , ..., d imply summationover them. Using the generators Γ a of the su( N c ) al-gebra, the gauge field can be written as a fundamentalrepresentation matrix A µ = A aµ Γ a . We consider the fol-lowing theories: • ‘3D’ or ‘3+1D’: in d = 3 spatial dimensions, suchthat S Y M [ A ] = S D Y M [ A ] . (2) • ‘2D’ or ‘2+1D’: in d = 2 spatial dimensions, suchthat S Y M [ A ] = S D Y M [ A ] . (3) • ‘2D+sc’ or ‘Glasma-like 2+1D’: originally in d = 3spatial dimensions where no field depends on thecoordinate x . This results in a 2+1D theory withthe classical action S Y M [ A ] = QL (cid:0) S D Y M [ A ] + S Dφ [ φ ] (cid:1) (4)with an adjoint scalar field φ a ≡ A a having theaction S Dφ [ φ ] = − (cid:90) d x ( D abµ φ b )( D µac φ c ) , (5)with summation over µ = 0 , , D abµ = δ ab ∂ µ − gf abc A cµ . Thelength in the x direction L drops out of the clas-sical dynamics. We have also factored out a con-stant momentum scale Q and rescaled the gaugecoupling and all fields as g Q / → g , A Q − / → A , . . . , such that the momentum dimensions of the ac-tion, the coupling constant and the fields match theones for 2+1D. This results in [ S Y M ] = [ S D Y M ] =[ S Dφ ] = 0, [ g ] = 1 / A ] = [ φ ] = 1 / Q used in therescaling is in principle arbitary, but we will usea value related to the conserved energy density, asdiscussed in more detail below.We will refer to this theory as Glasma-like 2+1D,since the Glasma occurring at the initial stages atultrarelativistic heavy ion collisions is also invariantin one spatial dimension (rapidity) [2, 6, 8]. Unlikethe expanding Glasma usually employed in mod-els of the initial collision dynamics, the 2+1D the-ory we consider here is defined in a nonexpandingcoordinate system, and our initial condition corre-sponds to a positive (but small) longitudinal pres-sure P L >
0. Also in the expanding Glasma theexpansion becomes less important at later times,as reflected e.g. in the fact that P L becomes pos-itive at τ (cid:38) /Q s . Thus our Glasma-like systemcould be interpreted as a nonexpanding model ofthis later Glasma state. B. Spectral and statistical correlation functions
Here we will give a brief overview of the spectral andstatistical correlation functions measured in this work.For a more comprehensive introduction and descriptionof the methods used we refer the reader to [40].The statistical correlation function is defined as theanticommutator of two field operators (cid:104) AA (cid:105) jk ( x, x (cid:48) ) = 12 d A (cid:88) a (cid:68)(cid:110) ˆ A aj ( x ) , ˆ A ak ( x (cid:48) ) (cid:111)(cid:69) (cid:104) EE (cid:105) jk ( x, x (cid:48) ) = 12 d A (cid:88) a (cid:68)(cid:110) ˆ E ja ( x ) , ˆ E ka ( x (cid:48) ) (cid:111)(cid:69) , (6)with spatial components j, k = 1 , . . . , d , the dimension ofthe adjoint representation d A = N − x ≡ ( t, x ). Because of E j ( x ) = ∂ t A j ( x ),the correlators (cid:104) EE (cid:105) and (cid:104) AA (cid:105) are related via timederivatives. In analogy, the statistical correlation func-tion for scalar components of the Glasma-like theory aredefined as (cid:104) φφ (cid:105) ( x, x (cid:48) ) = d A (cid:80) a (cid:68)(cid:110) ˆ φ a ( x ) , ˆ φ a ( x (cid:48) ) (cid:111)(cid:69) and (cid:104) ππ (cid:105) ( x, x (cid:48) ) = d A (cid:80) a (cid:104){ ˆ π a ( x ) , ˆ π a ( x (cid:48) ) }(cid:105) . The scalar fieldsare related to the fields of the original theory via φ ≡ A and π ≡ E . In the classical limit the statistical correla-tion functions are easy to measure, since the anticommu-tator of Heisenberg field operators reduces to a productof two classical fields.The spectral function, on the other hand, is defined asthe commutator of two field operators ρ jk ( x, x (cid:48) ) = id A (cid:88) a (cid:68)(cid:104) ˆ A aj ( x ) , ˆ A ak ( x (cid:48) ) (cid:105)(cid:69) ˙ ρ jk ( x, x (cid:48) ) = id A (cid:88) a (cid:68)(cid:104) ˆ E aj ( x ) , ˆ A ak ( x (cid:48) ) (cid:105)(cid:69) , (7)and analogously for the scalar components of the Glasma-like theory. We will mostly study its time derivative ˙ ρ ,which we will refer to as the “dotted spectral function”.In the classical limit the commutator corresponds to theDirac bracket, which generalizes the Poisson bracket forsystems with constraints (such as the Gauss’ law). There-fore, a direct measurement of the spectral function is a quite involved task. However, the spectral function isintimately related to the retarded propagator G Rjk ( t, t (cid:48) , p ) = θ ( t − t (cid:48) ) ρ jk ( t, t (cid:48) , p ) . (8)The latter can be extracted numerically by using ourlinear response framework [39, 40]. There we introducea small instantaneous source j for different momentummodes on top of the classical fields A ( x ), E ( x ). Thesource generates small response fields a ( x ), e ( x ). Thesefollow linearized equations of motion constructed to sat-isfy the Gauss law constraint exactly. The spectral func-tion ρ and the dotted spectral function ˙ ρ are computedas correlation functions of j with a and e , respectively.In (8) we implied a spatial Fourier transform with re-spect to x − x (cid:48) . Apart from working in momentum space,it is beneficial to compute the correlation functions in fre-quency space. For that, a Fourier transform with respectto the relative time ∆ t = t − t (cid:48) for fixed central time¯ t = ( t + t (cid:48) ) / (cid:104) EE (cid:105) ( t, ∆ t, p ) and ˙ ρ ( t, ∆ t, p ) with respect to ∆ t .In practice, as discussed in Refs. [40, 45], we approximatethis by a numerically more efficient transform where thelower limit t is kept fixed, and one Fourier transformswith respect to an upper limit t + ∆ t : (cid:104) EE (cid:105) ( t, ω, p )= 2 (cid:90) ∞ d∆ t cos( ω ∆ t ) (cid:104) E ( t + ∆ t/ E ( t − ∆ t/ (cid:105) ( p ) ≈ (cid:90) ∆ t max d∆ t cos( ω ∆ t ) h (∆ t ) (cid:104) E ( t + ∆ t ) E ( t ) (cid:105) ( p ) , (9)and analogously for ˙ ρ . We justify this approximation bythe observation that our systems at sufficiently late timesdepend on t only weakly within a relative time windowof the order of the inverse plasmon mass ∼ /ω pl , whichis the relevant timescale for ∆ t . If not stated otherwise,the resulting curves are made smoother by employingstandard signal processing techniques, as in Ref. [42]. Tothis end we have introduced in Eq. (9) a Hann windowfunction h (∆ t ) = 12 (cid:18) π ∆ t ∆ t max (cid:19) (10)and use zero padding, which implies evaluating Eq. (9) atmore intermediate frequencies than computed using thediscrete Fourier transform. As will be demonstrated inSec. III C, these methods do not change the forms of thecorrelation functions considerably, while reducing back-ground ringing that results from a finite time window inthe Fourier transform.In this work, we will not Fourier transform the corre-lation functions (cid:104) AA (cid:105) and ρ for 2+1D theories directly,because these functions turn out to oscillate around anon-zero value in the time domain ∆ t . In stead, we ob-tain the frequency space spectral function by Fourier-transforming the time derivative ˙ ρ ( t, ∆ t, p ), and subse-quently dividing by ω . Note that ρ ( t, ∆ t, p ) and ρ ( t, ω, p )are odd functions in ∆ t and ω respectively, and are re-lated by a Fourier sine transform in stead of the cosinetransform (9). Thus we have explicitly ρ ( t, ω =0 , p ) = 0.However, as we will see in Sec. III E, the oscillationsaround a non-zero value at large ∆ t lead to the actualspectral function not approaching this limit smoothlywhen ω → (cid:104) EE (cid:105) T = v jT v kT (cid:104) EE (cid:105) jk or (cid:104) EE (cid:105) L = v jL v kL (cid:104) EE (cid:105) jk , with v T =( p y , − p x ) T /p and v L = ( p x , p y ) T /p .The (single-particle) distribution function is an impor-tant quantity that we can use to characterize the dynam-ical state of the system. We will use the same definitionsas in our previous publications [36, 40, 45] f E ( t, p ) = (cid:104) EE (cid:105) T ( t, ∆ t =0 , p ) pf π ( t, p ) = (cid:104) ππ (cid:105) ( t, ∆ t =0 , p ) p . (11)The employed temporal gauge leaves room for gaugetransformations that only depend on the spatial coor-dinates. Since all the correlation functions discussed inthis subsection are not manifestly gauge-invariant observ-ables, we remove the residual gauge freedom for equaltime measurements by fixing to Coulomb-type gauge ∂ j A j = 0 (cid:12)(cid:12) t at the time t of the measurement, as oftenemployed in classical-statistical gauge simulations [25–30]. In this physical gauge, gauge fields are always trans-versely polarized while electric fields may include longi-tudinal terms. For unequal-time correlators we also im-pose the same condition at the initial time t for the mea-surement. During the time evolution, the system thengradually shifts away from the gauge condition so that itis not exactly satisfied at the second measurement time t + ∆ t . However, as argued in [40], this effect is expectedto be small for sufficiently short relative times ∆ t (cid:28) t ,which we also assumed in order to perform the Fouriertransform (9) in a finite ∆ t interval. C. Self-similarity and nonthermal fixed points
Many highly occupied systems exhibit non-thermalfixed points (also referred to as universal classical at-tractors). These have been seen in classical non-Abeliangauge theories [25, 27, 28, 36], scalar theories [46–50], lon-gitudinally expanding systems [30, 51, 52] and recentlyalso in ultra-cold atom experiments [53–55]. A character-istic feature of these attractors is that the system forgetsthe details of its initial conditions and can usually be un-derstood using a very simple set of scaling functions andpower laws, which tremendously simplifies the theoreticaldescription. In our previous work [36], we established that puregauge theories also in 2+1D exhibit self-similarity, as de-fined by f ( t, p ) = ( Qt ) α f s (( Qt ) β p ) , (12)and extracted the universal scaling exponents β = − / , α = 3 β . (13)We observed there that the scaling exponents in the2+1D systems can be understood using relatively sim-ple kinetic theory arguments, even if a kinetic theorydescription is not expected to work quantitatively. Thisresult was one of the main motivations for the presentpaper, since it hints that a quasiparticle description of2+1D theories at high momenta might be possible.To understand the physical meaning of Eqs. (12)and (13), we define the time-dependent hard scale Λ( t )as the momentum scale that contributes the most tothe perturbative estimate of the energy density ε ∼ (cid:82) d p p f ( t, p ). Then the scaling relation (12) and thescaling exponents implyΛ( t ) ∼ Q ( Qt ) − β , m D ( t ) ∼ Q ( Qt ) β , ε = const , (14)where we included the expected scaling of the soft scalerepresented by the Debye mass from Eq. (A3). Thus, thehard scale grows and the soft scale decreases with time,while energy density is conserved.Throughout this work all the measurements are per-formed at sufficiently large values of time t to be in theself-similar regime. D. Initial conditions
For the 2+1D systems, we use the same initial condi-tion as in [36]. We consider weakly coupled g /Q (cid:28) f (cid:29) f ( t, p ) at the initial time t = 0 both for gauge and scalar excitations (i.e., f E = f π = f ) is given by f ( t = 0 , p ) = Qg n e − p p . (15)Here Q is a gauge invariant momentum scale, which ismore precisely defined by Q ≡ (cid:115) C g εd pol ( N − . (16)We consider systems in a fixed size box, and thus theenergy density is a conserved quantity. Since on theclassical level it is possible to carry out field simulationswithout making explicit reference to the precise value ofthe coupling g , the quantities we have access to are theenergy density scaled with the coupling g ε , which hasthe momentum dimension [ g ε ] = 4, and g f . Unlessstated otherwise, we will use the initial occupation num-ber n = 0 . , for which p = Q for our chosen definitionas detailed below. We have previously shown [36] thatthe form of the initial conditions is unimportant since thesystems will approach an attractor solution that only de-pends on Q in the sense of Eq. (12). The number of non-longitudinal polarizations in (15) is d pol , with d pol = 1 forthe 2+1D and d pol = 2 for the Glasma-like 2+1D theo-ries. The constant C is taken as C = 20 √ π ≈ , whichwe merely choose for convenience such that for n = 0 . p = Q .In 2+1D the coupling constant g is dimensionful: if onekeeps the dimensionless combination g /Q constant, oneobserves that (16) leads to the proportionality Q ∝ √ ε, which is natural for a scale derived from a 2-dimensionalenergy density. Combining the definition (16) with a per-turbative estimate ε ≈ d pol ( N − (cid:82) d p/ (2 π ) p f ( t =0 , p ) for the energy density, one obtains Q ≈ C (cid:90) d p (2 π ) p g f ( t = 0 , p ) Q = 10 n p , (17)which is independent of g , d pol and N c .For the 3+1D theory we employ the same isotropicinitial condition as in [40] with the distribution functioninitially given by f ( t = 0 , p ) = n g p p e − p p . (18)Here the coupling is dimensionless. In this work we use n = 0 . Q = √ n p ∝ (cid:112) g ε .In our figures, all dimensionful quantities are rescaledby appropriate powers of Q of the respective theory tomake them dimensionless, unless a different rescaling pre-scription is stated explicitly. E. HTL expressions
Diagrammatically, the HTL approximation corre-sponds to an all-order resummation with the kinematicapproximation that the external lines are considered tobe soft compared to the hard momenta flowing in theinternal lines of the diagrams. This corresponds to thenon-Abelian generalization of the Vlasov equations, theWong equations, where the soft modes are consideredto be classical fields and the hard modes correspond toclassical particles [38]. While we expect the interactionbetween the soft and hard modes to be described bythe HTL theory as is the case in 3+1D, the paramet-ric counting for 2+1D in [36] suggests that in additionto the interaction between soft and hard modes, there is also a leading-order interaction between the soft modesamong themselves, which is not captured by the HTL ap-proximation. Therefore, we expect deviations from HTLexpressions arising from the nonperturbative soft-soft in-teractions. We compare our numerical results to the HTLcalculations to quantify the deviations.Here we summarize the main HTL results that are rel-evant for the comparison with our data. Details are writ-ten in Appendix A. The HTL spectral function can bedecomposed into a transversely polarized, a longitudi-nally polarized and, for the Glasma-like theory, a scalarcontribution, denoted below by the index α = T , L , or φ ,respectively. Each contribution can be further split intoa Landau damping part for | ω | < p and a quasiparticlepart ρ HTL α ( ω, p ) = ρ Landau α ( ω, p )+ 2 πZ α ( p ) (cid:2) δ (cid:0) ω − ω HTL α ( p ) (cid:1) + δ (cid:0) ω + ω HTL α ( p ) (cid:1)(cid:3) . (19)The Landau damping expressions in 2+1D read ρ Landau T = 2 m D x √ − x θ (1 − x )((1 − x ) ( p/m D ) + x ) + x (1 − x )(20) ρ Landau L = 2 x m D / √ − x θ (1 − x )( p + m D ) + x m D / (1 − x ) (21) ρ Landau φ = 0 , (22)with x = ω/p . The Debye mass m D entering the ex-pressions is determined within HTL at leading orderin Eq. (A3) and depends on the distribution function f ( t, p ). It is connected with the asymptotic mass via m D = m in 2+1D theories and m D = 2 m for the3+1D system. In this work, we compute the asymptoticmass like in [40] in d spatial dimensions as m = d pol N c (cid:90) d d p (2 π ) d g f ( t, p ) (cid:112) m + p , (23)which is a self-consistent generalization of Eq. (A3) andreduces the dependence on the definition of the distribu-tion function [40].The dispersion relations in 2+1D are ω HTL T ( p ) = m HTL (cid:115) p − p + (cid:112) p − p (24) ω HTL L ( p ) = m HTL p (cid:112) p (25) ω HTL φ ( p ) = (cid:113) m + p , (26)where we defined ˜ p = p/m HTL . The expressions for thequasiparticle residues Z α ( p ) are written in Appendix A. III. NUMERICAL RESULTS
In this section, we present our numerical results forthe excitation spectrum of non-equilibrium 2+1D the-
FIG. 1. The extracted and normalized transverse sta-tistical correlation function (cid:104) EE (cid:105) T ( t, ω, p ) / (cid:104) EE (cid:105) T ( t, ∆ t =0 , p )is shown for a) the 2 + 1 dimensional theory, b) for theGlasma-like theory and c) for an isotropic 3 + 1 dimen-sional system. Additionally, we show the scalar correlator (cid:104) ππ (cid:105) ( t, ω, p ) / (cid:104) ππ (cid:105) ( t, ∆ t =0 , p ) of the Glasma-like system in d).Note that all amplitudes have been cut off at the same fixedvalue 22 /Q corresponding to the red region. For compari-son, we added the HTL dispersion relations ω HTL
T/φ ( p ) as blackdashed lines and a relativistic dispersion ω rel ( p ) = (cid:113) ω + p as a continuous gray line with the extracted values for ω pl asin the text. ories by extracting spectral and statistical correlationfunctions. We briefly discuss some results in (classical)thermal equilibrium in Appendix B.For our simulations, we use the lattice size 1024 andlattice spacing Qa s = 1 / Qa s = 1 / lattice for the Glasma-like 2+1Dsystem. We showed in [36] for equal-time correlationfunctions that for the considered initial conditions, thesediscretization parameters are sufficient to avoid latticeartifacts. Our resulting correlation functions will be com-pared with those from isotropic 3+1D systems computedin [40] with Qa s = 0 . lattice, which has beendemonstrated not to show any significant lattice artifactsfor this choice of parameters.While respecting the Gauss law constraint to almostmachine precision, the linearized fluctuations on top ofthe inhomogenous, time-dependent, background field donot have conservation laws of quantities like energy, mo-mentum and angular momentum. Thus, there is no guar-antee that they would not grow strongly with time, lead-ing to increasing error bands. Indeed, we observe sucha growth, more strongly in two than in three spatial di-mensions. This limits in practice our ability to extract FIG. 2. The normalized longitudinal statistical correlator (cid:104) EE (cid:105) L ( t, ω, p ) / (cid:104) EE (cid:105) L ( t, ∆ t =0 , p ) is plotted for a) the 2 + 1dimensional theory, b) for the Glasma-like theory and c) foran isotropic 3 + 1 dimensional system. We also included theHTL dispersion relation ω HTL L ( p ) as black dashed lines, a rel-ativistic dispersion ω rel ( p ) as a continuous gray curve and thefree dispersion ω = p as gray dash-dotted lines. the spectral function in the time domain. In practicewe do not extend our extraction of the spectral functionfor 2+1D theories to more than ∆ t ≤ ∆ t max = 80 /Q .No such restrictions apply for statistical correlation func-tions, which allows us to simulate to much later times.While it suffices to use ∆ t max = 80 /Q also for transverseand longitudinal statistical correlators in 2+1D theories,unless stated otherwise, we employ ∆ t max = 200 /Q forthe scalar statistical correlator of the Glasma-like theoryto correctly capture the long-lived scalar excitations atlow momenta.All the results in the figures are shown at Qt = 500 forthe 2+1D systems and Qt = 1500 for the 3+1D theory,unless stated otherwise. From our earlier papers [36, 40]we see that the systems at these times are well in theuniversal scaling regime discussed in Sec. II C. We willdiscuss how the correlators depend on time t in Sec. III F.Correlation functions as functions of frequency or rel-ative time for fixed momenta are shown with error bars.These result from statistical averaging of the correlationfunctions over 320 - 400 configurations in frequency spaceor in the relative time domain, respectively. A. Transverse, longitudinal and scalar correlations
Figure 1 shows the extracted normalized statistical cor-relation functions (cid:104) EE (cid:105) α ( t, ω, p ) / (cid:104) EE (cid:105) α ( t, ∆ t =0 , p ) asfunctions of frequency and momentum for transverse ex-citations in 2+1D theory (a), Glasma-like 2+1D (b),and isotropic 3+1D (c), and for scalar excitation in theGlasma-like 2+1D theory (d).All frequencies and momenta of the correlators are nor-malized using the mass ω pl , which is the plasmon fre-quency of gluonic polarizations. While there are multi-ple ways of determining ω pl [56, 57], we extract it fromthe dotted spectral function ˙ ρ T ( t, ω, p =0) at vanishingmomentum using a (Gaussian) fit function, as explainedin Sec. III C. At the considered times, the values of theplasmon frequency read ω D pl = 0 . Q for the 2+1D the-ory, ω D + sc pl = 0 . Q for the Glasma-like 2+1D theoryand ω D pl = 0 . Q for the 3+1D theory. The continuousgray lines in Fig. 1 correspond to a relativistic disper-sion ω rel = (cid:113) ω + p with the numerically extracted ω pl . We show the transverse and scalar HTL dispersionrelations ω HTL
T/φ ( p ) from Eq. (24) and Eq. (26) as blackdashed lines. Note that the relativistic dispersion is anad hoc ansatz with the parameter ω pl extracted by fittingthe data. The HTL dispersion relation, in contrast, has afunctional form determined by the leading order HTL cal-culation, and a mass scale determined by the distributionfunction f ( t, p ) using Eq. (23). For the HTL mass, theresulting values are m D HTL = 0 . Q , m D + sc HTL = 0 . Q and m D HTL = 0 . Q for the considered parameters andtimes of the different theories.When considering the correlation functions in Fig. 1for fixed momentum, one observes an excitation peakas a function of frequency for each of the different non-Abelian gauge theories and their transverse gluonic orscalar contributions. The additionally shown relativisticdispersion ω rel ( p ) and HTL dispersion relations lie wellwithin the width of the peak. Interestingly, the scalar ex-citation shows some deviations from a relativistic disper-sion, which actually is the leading order HTL prediction(26). These deviations are at their strongest at p ∼ ω pl .However, at larger momenta p (cid:29) ω pl the dispersion rela-tion coincides with its HTL prediction, which also agreeswith the transverse HTL dispersion for large momenta.As we will see in Sec. III D, also the shape of the spectralfunction agrees between transverse and scalar polariza-tions at higher momenta.An important difference to the expected form of theexcitations within the hard-loop framework lies in thewidth of the peaks (considered as functions of ω ), whichin hard loop theory is a next-to-leading order effect andis thus supposed to be small. We also observe that forgluonic excitations the peaks in 2+1D theories are muchwider than in 3+1D or for low-momentum scalar excita-tions, whose width seems to be comparable to the gluonicwidth in 3+1D. Transverse gluonic excitations of 2+1D Note that in the lower panels of Fig. 1, the amplitudes of the3+1D and scalar excitations have been cut off at the same value theories at low momenta p (cid:46) ω pl also show sizeable con-tributions all the way to ω = 0. We will discuss thisobservation in more detail below in Sec. III E.Figure 2 shows the normalized longitudinal statisti-cal correlator (cid:104) EE (cid:105) L ( t, ω, p ) / (cid:104) EE (cid:105) L ( t, ∆ t =0 , p ) of 2+1D(a), Glasma-like 2+1D (b), and isotropic 3+1D theory (c)for comparison. The gray dash-dotted lines correspondto the free dispersion ω = p , the gray continuous lines to ω rel ( p ) and the black dashed curves to the HTL disper-sion relations ω HTL L ( p ) from Eq. (25). One observes thatthe peaks at low momenta p (cid:46) ω pl have a similar widthand magnitude as the corresponding transversely polar-ized gluonic peaks of Fig. 1, and agree with the longitu-dinal HTL dispersion reasonably well. At all momenta,one finds that the longitudinal correlations involve a fi-nite valued continuum for ω (cid:46) p , which corresponds tothe Landau cut contribution. Different from transversedotted spectral functions, the quasiparticle peaks of lon-gitudinal correlators are strongly suppressed at highermomenta p (cid:38) ω pl and Landau damping becomes thedominant contribution, as visible in the figure.These properties will be studied in more detail by com-paring with HTL calculated curves in the following sub-section. Note that while we have shown the statisticalcorrelation functions (cid:104) EE (cid:105) α ( t, ω, p ) here, they appear tobe related to the respective dotted spectral functions˙ ρ α ( t, ω, p ) by the generalized fluctuation-dissipation re-lation (27) even in our far-from-equilibrium situation, aswill also be shown in the following. Therefore, the dis-cussions of this subsection are also valid for the dottedspectral functions in the considered theories. B. Comparison of spectral and statisticalcorrelation functions
In the isotropic 3+1D theory, we observed in [40] thatthe spectral and statistical correlation functions in ourconsidered far-from-equilibrium system are related by thegeneralized fluctuation-dissipation relation (cid:104) EE (cid:105) α ( t, ω, p ) (cid:104) EE (cid:105) α ( t, ∆ t =0 , p ) ≈ ˙ ρ α ( t, ω, p )˙ ρ α ( t, ∆ t =0 , p ) . (27)We will show here that this relation also holds in 2+1Dtheories for each polarization α = T, L, φ . The respec-tive equal-time dotted spectral functions are determinedby the sum rules (A15)-(A17) while the equal-time sta-tistical correlators (cid:104) EE (cid:105) α ( t, ∆ t =0 , p ) are extracted nu-merically. The generalized fluctuation-dissipation re-lation (27) implies that the normalized statistical anddotted spectral functions agree quantitatively, which as for the transverse 2+1D excitations in the upper panels toenable a direct comparison. However, their amplitudes reachvalues around 80, such that the peak width visible in the figureappears larger than in reality. p = 0.049 Q ρ . T
FIG. 3. The transverse normalized statistical correlator (cid:104) EE (cid:105) T ( t, ω, p ) / (cid:104) EE (cid:105) T ( t, ∆ t =0 , p ) and dotted spectral func-tion ˙ ρ T ( t, ω, p ) for both 2+1-dimensional theories at differentmomenta as functions of frequency. They lie on top of eachother within numerical precision, demonstrating the validityof the generalized fluctuation-dissipation relation (27). TheHTL expressions for ωρ are shown as green dashed lines whileblack lines denote fits using the Gaussian distribution in (30). is nontrivial out of equilibrium. In thermal equilib-rium both correlations are connected by the fluctuation-dissipation relation, which for low momenta p (cid:28) T be-low the temperature reads (cid:104) EE (cid:105) α ( ω, p ) = T ˙ ρ α ( ω, p ).In contrast, out of equilibrium their connection may bemuch more complicated or even absent, as observed forscalar theories [41, 42, 58]. The HTL framework pre-dicts a generalized fluctuation-dissipation relation [59]with a time-dependent effective temperature T eff ( t ) = (cid:104) EE (cid:105) α ( t, ∆ t =0 , p ) / ˙ ρ α ( t, ∆ t =0 , p ) (cid:104) EE (cid:105) α ( t, ω, p ) ≈ T eff ( t ) ˙ ρ α ( t, ω, p ) (28)to hold for sufficiently small momenta p (cid:28) Λ( t ). HereΛ( t ) is a time-dependent hard scale that dominates theenergy density and plays the role of a temperature scalefor momenta in the non-equilibrium setting. However,we emphasize that we observe the relation (27) to besatisfied in general. This also includes higher momenta p (cid:38) Λ( t ), where (cid:104) EE (cid:105) α ( t, ∆ t =0 , p ) / ˙ ρ α ( t, ∆ t =0 , p ) arenon-constant functions of t and p .The validity of the generalized fluctuation-dissipationrelation (27) in 2+1D non-Abelian plasmas is one of themain results of our work. It is demonstrated in Figs. 3and 4, which show the normalized statistical and spec-tral correlation functions for three different momenta forthe 2+1D (left panels) and for the Glasma-like theory p = 0.049 Q ρ . L
FIG. 4. The longitudinal normalized statistical corre-lator (cid:104) EE (cid:105) L ( t, ω, p ) / (cid:104) EE (cid:105) L ( t, ∆ t =0 , p ) and normalized dot-ted spectral function ˙ ρ L ( t, ω, p ) / ˙ ρ L ( t, ∆ t =0 , p ) for both 2+1Dtheories for the same momenta and with analogous HTL andfitting curves as in Fig. 3. For p = 0 . Q we used the smallertime window ∆ t max = 60 /Q to reduce ringing of the dottedspectral function. (right panels) for transverse and longitudinal polariza-tions, respectively. The normalized statistical correla-tion functions (blue points) precisely overlap with thecorresponding dotted spectral functions (red points), upto residual ringing of the dotted spectral functions thatresults from the finite window size and statistical fluctua-tions of the Fourier transform. Thus, the figures confirmthe generalized fluctuation-dissipation relation given byEq. (27).We also observe in these figures that the width andshapes of the correlation functions look very similar inboth the 2+1D and Glasma-like 2+1D theories for trans-verse and longitudinal excitations, when plotted as func-tions of ω/Q . Note that the scale Q is related to theenergy density per polarization state . Thus, adding thescalar degrees of freedom to the system without chang-ing the distribution f ( p ) of the gauge particles does notseem to modify the shape or the width. In contrast tothe width of the peak, the mass in the Glasma-like sys-tem is always larger, as seen particularly well by the shiftof the peak at lower momenta. Thus, the squared plas-mon mass is proportional to the number of polarizationstates, as expected from HTL in Eq. (23). Taking to-gether these two observations seems to indicate that, toa first approximation, the scalar sector contributes to themass but does not contribute to the damping of gluonicexcitations.If we plotted the correlation functions at the same p/ω pl as functions of ω/ω pl , i.e., normalized by theplasmon mass of the corresponding theory instead ofan energy-related momentum scale, the peak positionswould overlap by construction. However, the peaks inthe 2+1D Glasma-like system would appear more narrowthan in the 2+1D theory. This effect, which we have al-ready seen in Figs. 1 and 2, means that in the Glasma-liketheory the quasiparticle excitations are slightly longerlived than in pure 2+1D.We now compare the correlators to the transverse andlongitudinal spectral functions calculated in HTL (19),which consist of a Landau cut contribution and a quasi-particle delta peak. They are shown as green dashedlines in Figs. 3 and 4. Importantly, different from 3+1Dgluonic plasmas, the extracted dotted spectral functionsare so broad that an accurate distinction between Lan-dau cut region and an excitation peak becomes difficultfor both polarizations. Therefore, the HTL expressionsmostly do not provide a good description of the nonper-turbative simulation results of the 2+1D systems.One exception is the large-momentum region p (cid:38) m HTL of the longitudinally polarized dotted spectralfunction, which concerns the central and lower panels ofFig. 4. According to HTL expectations that are discussedin Appendix A, the Landau damping term is expected todominate the spectral function for these momenta. Weindeed observe that the HTL Landau damping curves forlow frequencies ω (cid:28) p agree well with the numerical data.This observation, which validates the HTL formulas forthis specific case, is nontrivial, since the HTL frameworkis not expected to work well for the 2+1D theories dueto the importance of the missing soft-soft interactions, asdiscussed in Sec. II E.Note that the sharp peaks of the HTL Landau con-tributions around ω ∼ p are smoother in the numeri-cally extracted correlators. We have reported of a similarsmoothening of the Landau cut region for 3+1D plasmasin [40], which may result from including more mode in-teractions than considered in HTL at leading order. C. Shape of the peaks
Let us now discuss the shape of the excitations inmore detail. Computations within the HTL formalism(Sec. II E) suggest that the normalized correlation func-tions can be written as a sum˙ ρ α ( t, ω, p )˙ ρ α ( t, ∆ t =0 , p ) ≈ h ( ω, p ) + ω ρ Landau α ( t, ω, p )˙ ρ α ( t, ∆ t =0 , p ) (29)with an excitation peak h ( ω, p ) that originates from apole at frequency ± ω α ( p ) with residue Z α ( p ) and be-comes a Delta function at leading order in HTL. A fi-nite peak width should thus be a subleading-order effect.Since the peak results from poles in the retarded prop-agator (Eqs. (A10)-(A12)), its shape could be expectedto be of Breit-Wigner form, which becomes a Lorentzian ω / Q N o r m . c o rr e l a t i on < EE > T Q Δ t max = 150Q Δ t max = 80no windowing 0 0.2 0.4 0.6Q Δ t max = 150Gauss fi tSech fi tBW fi t FIG. 5. The normalized transverse statistical correlator (cid:104) EE (cid:105) T ( t, ω, p ) / (cid:104) EE (cid:105) T ( t, ∆ t =0 , p ) of Glasma-like 2+1D sys-tem for momentum p = 0 . Q . In the left panel, differentways to compute the Fourier transform are compared, whichare seen to agree well. In the right panel, we compare differ-ent fit functions for the peak shape, as explained in the maintext. distribution for narrow peaks. A Lorentzian shape forthe excitation peaks has been found in 3+1D gluonic sys-tems [40] and would usually correspond to quasiparticleexcitations.The peak form is investigated in Fig. 5. In the rightpanel, the normalized transverse statistical correlator forthe larger time window Q ∆ t max = 150 is shown in redfor the Glasma-like theory at momentum p = 0 . Q . It iscompared to fits h ( ω, p ) that are modeled as a Gaussian,hyperbolic secant and Breit-Wigner distribution, respec-tively, given by h Gauss = Aγ (cid:114) π (cid:34) − (cid:18) ω − ω R γ (cid:19) (cid:35) + [ ω R (cid:55)→ − ω R ](30) h Sech = Aγ π (cid:20) π ω − ω R γ (cid:21) + [ ω R (cid:55)→ − ω R ] (31) h BW = a ω ( ω − ω R − γ ) + 4 ω γ (32)with sech( x ) = 1 / cosh( x ). Here [ ω R (cid:55)→ − ω R ] denotesthe previous term with the indicated substitution. Weneglect the Landau cut contribution in the fitting pro-cedure because it is subleading and does not describeour data well at frequencies ω ∼ p , which lie within thewidth of the excitation peak. Comparing the differentfits in the plot, one observes that the Breit-Wigner dis-tribution fails to describe the peak faithfully, which isparticularly well visible at the right tail. In contrast, theshape of the peak is well described by the non-Lorentziandistributions h Gauss and h Sech . Indeed, we have included All parameters, which are the dispersion relations ω D and thewidth γ are momentum and time dependent and will be discussedin Sec. III F. The amplitude A corresponds to the residue of thequasiparticles and is close to unity for transverse excitations. h Gauss for different momenta and polarizations alreadyin the Figs. 3 and 4 as black lines. The fits are seen toagree well with the shape of the dotted spectral functionsfor all momenta, where such a comparison is sensible.Note that for longitudinal spectral functions, quasiparti-cle excitations are expected to be suppressed at highermomenta p (cid:38) m D . Therefore, we only included fits forsmall momenta p (cid:28) m D in Fig. 4.We have checked that the non-Lorentzian shape is notan artefact of a finite time window in the Fourier trans-form or of signal processing techniques. In the left panelof Fig. 5 we show the same correlation function as inthe right panel for the time windows Q ∆ t max = 150and Q ∆ t max = 80 with Hann windowing as well asfor Q ∆ t max = 80 without windowing, which implies h (∆ t ) = 1 in Eq. (9). One observes that all curves accu-rately agree within uncertainties.This non-Lorentzian peak, thus, appears to be a dis-tinct feature of 2+1D theories in contrast to 3+1Dplasmas, where Lorentzian excitations emerge instead[40]. A similar non-Lorentzian shape has been encoun-tered in single-component non-relativistic and relativis-tic scalar models at low momenta [41, 42]. It was alsoused to distinguish the corresponding excitations fromthe usual Lorentzian peaks that dominate in O ( N )-symmetric scalar models for a large number of compo-nents N (cid:29) D. Scalar excitation
We can also compare correlation functions in the (rel-ative) time domain. Consistently with Eq. (27), the gen-eralized fluctuation-dissipation relation (cid:104) EE (cid:105) α ( t, ∆ t, p ) (cid:104) EE (cid:105) α ( t, ∆ t =0 , p ) ≈ ˙ ρ α ( t, ∆ t, p )˙ ρ α ( t, ∆ t =0 , p ) , (33)is also satisfied in ∆ t , where α denotes the polariza-tions T, L, φ . This is demonstrated in the upper pan-els of Fig. 6 for the 2+1D Glasma-like theory, where thenormalized statistical (cid:104) EE (cid:105) α ( t, ∆ t, p ) / (cid:104) EE (cid:105) α ( t, ∆ t =0 , p )and (dotted) spectral correlation functions ˙ ρ α are seento nicely coincide for scalar and transverse polarizationsseparately. Interestingly, both polarizations even agreewith each other at high momenta p (cid:29) m D (right panels),while they are seen to differ at lower momenta p (cid:46) m D (left panels).This observation implies that the distinction betweentransverse and scalar excitations is only relevant for lowmomenta. We confirm this interpretation in the fre-quency domain in the lower panels of Fig. 6. In theright panel we show the correlation functions for thesame higher momentum and also include a (Gaussian)fit to the scalar correlator given by (30). All correla-tors lie on top of each other, revealing the same shape, -1-0.500.51 0 50 100 150 Time: Q Δ tp = 0 Q N o r m . c o rr e l a t i on s ππ > ρ . ϕ
The correlation functions (cid:104) ππ (cid:105) / (cid:104) ππ (cid:105) (∆ t =0 , p ), ˙ ρ φ , (cid:104) EE (cid:105) T / (cid:104) EE (cid:105) T (∆ t =0 , p ) and ˙ ρ T asfunctions of time ∆ t . Note that in the left panel the rangeof ∆ t is longer. Bottom:
The same correlators as in the toprow but now in the frequency domain ω . The scalar dottedspectral function ˙ ρ φ is not shown because Q ∆ t max of at least200 would be necessary while the correlator becomes verynoisy for (cid:38)
80. Gaussian fits to the scalar correlator (cid:104) ππ (cid:105) areshown as black dashed lines. dispersion and width among transverse and scalar polar-izations. In contrast, at the lowest momentum p = 0, aclear difference is visible between the gluonic (cid:104) EE (cid:105) T andthe scalar correlation (cid:104) ππ (cid:105) . The scalar excitation is muchmore narrow, leading to longer-lived quasiparticles thanfor transverse gluonic excitations. We have observed thisproperty also in Fig. 1d) for low momenta p (cid:46) m D , whichis similar to the gluonic excitations in a spatially three-dimensional system. Furthermore, the peak position isshifted towards a slightly higher frequency. Interestingly,this agrees qualitatively with the HTL dispersion rela-tions (Appendix A 3), which predict that the peak posi-tions of the gluonic and scalar excitations at p = 0 are ω pl and m D , respectively, with ω pl < m D . E. Low frequency behavior
We now turn to the actual spectral functions ρ α , onlyhaving discussed their time derivative ˙ ρ α so far. For thescalar polarization, we see in Fig. 6 that ˙ ρ φ as a functionof ∆ t always oscillates around zero and has a narrow1 -6-4-20246 0 20 40 60Time: Q Δ tp = 0 Q S pe c t r a l f un c t i on s ρ ϕ ρ T FIG. 7. The actual scalar and transverse spectral functions ρ φ and ρ T of the Glasma-like 2+1D theory for momenta p =0 Q (left) and p = 0 . Q (right) as functions of relative time∆ t . Different from the dotted spectral function ˙ ρ T in Fig. 6, ρ T is seen to oscillate around a non-zero value in the leftpanel. p = 0 Q Spectral function ρ T p = 0.1 Q Spectral function ρ T p = 0 Q Frequency: ω / Q Spectral function ρ T Glasma-like 2+1D0 0.2 p = 0.1 Q
Frequency: ω / Q Spectral function ρ T FIG. 8. Transverse spectral function ρ T computed as ˙ ρ T /ω for small momenta p = 0 Q and p = 0 . Q as a function of ω for2+1D (left) and Glasma-like 2+1D (right) theories. Due toits symmetry, the spectral function satisfies ρ T ( t, ω =0) = 0identically while it follows ρ T ∼ /ω for ω → ρ T . excitation peak in the frequency domain with a vanish-ing value for ω →
0. As visible in Fig. 7, the scalarspectral function ρ φ also oscillates around zero. Thisconfirms that its Fourier transform ρ φ ( t, ω, p ) smoothlyapproaches zero as ω → ρ T ( t, ω, p ) is seen to behave differently thanthe scalar one. This is visible in frequency domain inFigs. 3 and 6, where the spectral function approachesnon-vanishing values for ω → p (cid:46) m D .As a consequence, the resulting spectral function shownin Fig. 8 (computed as ρ T = ˙ ρ T /ω as discussed ear-lier) behaves as ∼ /ω for ω →
0. Thus, although ρ T ( t, ω =0) = 0 due to the odd symmetry, the spectralfunction actually behaves as ρ T ( t, ω ) ∼ /ω for ω → ω . This is consistent withthe behavior of ρ T as a function of ∆ t for low momenta. This can be seen in the left panel of Fig. 7 where we ob-serve that ρ T does not oscillate around zero but possiblyapproaches a nonzero value for ∆ t → ∞ . This behav-ior is surprising since, based on the analytical structureof the HTL self-energies, we would expect ˙ ρ T and ρ T tosmoothly vanish as ω → ρ T ( t, ∆ t, p )to oscillate around zero.Such a behavior could, for instance, follow from havingan additional excitation at ω ≈ ρ T ( t, ω =0 , p ) there for low momenta. Thus, the smallviolation of time reversal invariance in our system doesnot seem to be the explanation either.Another possibility is that this observed feature is asso-ciated with our gauge fixing procedure that is discussedin the end of Sec. II B. Assuming ∆ t (cid:28) t , we fix to aCoulomb-type gauge at the time t , but the system doesnot exactly remain in Coulomb gauge at t + ∆ t . Tocheck whether this approximation may lead to such aneffect, one could perform the calculation fully in Coulombgauge, which would avoid the mentioned approximationbut require reintroducing a temporal component for thegauge potential. Moreover, it would be interesting to lookat the heavy quark diffusion coefficient (as in Ref. [45]),which is a gauge invariant observable sensitive to thesame physical scales. We leave these studies to a furtherpaper. To summarize, we do not have a clear interpreta-tion of the finite value of ˙ ρ at ω = 0. F. Time dependence of dispersion relations &damping rates
Now we will study how the dispersion relations ω α ( t, p )and excitation widths γ α ( t, p ) of the 2+1D systems de-pend on the time t . We recall that in the scaling solution(see Sec. II C) this can be used to obtain informationabout the dependence on the soft and hard scales, whichfollow separate power laws in t . Figure 9 shows our nu-merical results for the dispersion relations and dampingrates at different times Qt for the transverse gluon ex-citation in the 2+1D theory and for the transverse andscalar excitations of the Glasma-like 2+1D theory. Theywere extracted by fitting h Gauss ( ω, p ) to the normalizedtransverse gluon and scalar excitations, as explained inSec. III C. The main observation here is that all curveslie on top of each other when plotted in terms of thetime dependent mass ω pl ( t ). For the dispersion relationthis is relatively trivial, and just shows that the func-tional form of the dispersion relation does not change intime and only depends on the current value of the mass.We have added HTL and relativistic dispersion relationsas black dashed and continuous gray lines, respectively.2 ω T / ω p l ( t ) Qt = 200Qt = 500Qt = 2000 00.20.40.60.81 0 2 4 6 γ T / ω p l ( t ) γ T / Q Glasma-like 2+1D, transverse ω T / ω p l ( t ) Qt = 200Qt = 500Qt = 2000 00.20.40.6 0 2 4 6 γ T / ω p l ( t ) γ T / Q ω pl (t) Glasma-like 2+1D, scalar ω ϕ / ω p l ( t ) Qt = 200Qt = 500Qt = 2000 00.20.40.6 0 2 4 6 γ ϕ / ω p l ( t ) γ ϕ / Q FIG. 9. The peak position (dispersion relation) ω α ( t, p ) andthe peak width (damping rate) γ α ( t, p ) as functions of momen-tum p at different times for ( top :) the transversely polarizedgluons of the usual 2 + 1-dimensional theory and for ( center :)the transversely polarized gluons and ( bottom :) the scalarcorrelator of the Glasma-like 2 + 1-dimensional theory. Dis-persion relations, widths and momenta are normalized by thetime dependent mass scale ω pl ( t ) in the main plots and by thetime independent energy scale Q in the insets. Additionally,we show, in the left panels, the same HTL dispersion rela-tions ω HTL
T/φ ( p ) and relativistic dispersions ω rel ( p ) as in Fig. 1as black dashed and continuous gray lines, respectively. We observe that both agree sufficiently well with the fit-ted values to the numerical data. Small deviations liewithin the width of the excitation peaks, as discussed inSec. III A and are visible in Fig. 1.The Qt depedence of the width, shown in the rightpanels of Fig. 9, is much more interesting. Although thescaling is not quite as good as for the dispersion relation,the plots show that both the p -dependence and the mag-nitude of the damping rate follow the time dependence ofthe plasmon mass. This is different from the 3+1D case,where γ ( t ) was found to decrease faster than ω pl ( t ) withtime [40]. Thus, while for the 3+1D theory the quasipar-ticle peaks get narrower as the scale separation between ω / ω pl (t) ω p l ( t ) ρ . ( t, ω , p = ) Qt = 200Qt = 500Qt = 2000 0 1 2 3 4
Glasma-like 2+1D, p=0
Qt = 200Qt = 500Qt = 2000
FIG. 10. The dotted spectral function of gluonic excita-tions as a function of frequency for vanishing momentum ω pl ˙ ρ ( t, ω/ω pl , p =0) at different times, rescaled by the plas-mon frequency such that they are made dimensionless. Allcurves fall on top of each other, showing the correlation func-tions follow a self-similar evolution with the only parameter ω pl ( t ). the dynamical hard and soft scales grows with increasingtime, this does not happen in the 2+1D theories: thewidth of the excitation peak relative to the mass remainsconstant.As a consequence of this scaling, correlators rescaledwith appropriate powers of ω pl ( t ) remain time-independent as functions of frequency ω/ω pl ( t ) and mo-mentum p/ω pl ( t ). This is demonstrated in Fig. 10 at theexample of ˙ ρ T at p = 0 for both 2+1D systems. Overa wide range of times Qt = 200 , , ω pl ( t ). In fact, also in 2D classical thermalequilibrium the width of the excitations is of the order ofthe mass, as we show in Appendix B. This is suggestiveof a general relation γ α ∼ ω pl for 2+1D systems.The behavior γ α ∼ ω pl should be contrasted with theexpectation from HTL, where the width is expected tobe proportional to the “effective temperature of the softmodes” g T ∗ ( t ). Indeed, in the 3+1D case we haveobserved [40] that the time dependence of γ α is con-sistent with this expectation. For the 2+1D theoriesa simple parametric estimate would lead to the scaling g T ∗ ∼ Q ( Qt ) − / [36], which is a much stronger time de-pendence than what we observe since ω pl ∼ Q ( Qt ) − / .Thus, the behavior of γ α ( t ) in the 2+1D theories indi-cates that the phenomenon of plasmon decay happens ina qualitatively different way than in 3+1D.As a consequence, excitations at low momenta p (cid:46) m D ( t ) for different times remain too short-lived to beconsidered as quasiparticles. Therefore, while an effectivekinetic theory description may exist for larger momenta,its collision kernel, depending on the soft modes, mustbe determined nonperturbatively.3 IV. CONCLUSION
In this paper we have studied the excitation spectrumof 2+1D classical gluodynamics in the self-similar regime.Our aim was to understand whether the boost invariantsystem created at the initial stages of ultrarelativisticheavy-ion collisions can be understood in a quasiparticlepicture. We have studied two theories. The first one is agenuine 2+1D theory where the longitudinal direction iscompletely absent. The second theory is Glasma-like, inthe sense that it is obtained from a 3+1D theory that isinvariant in the longitudinal direction. As a consequence,the gauge field in this direction becomes a scalar field.One of our main observations is that excitations arebroader in 2+1D than in isotropic 3+1D systems. Thisis seen in the damping rates, which correspond to thewidths of the excitation peaks in statistical and spectralcorrelation functions. In the 3+1D theory, the dampingrate is a subleading effect, and an increasing scale separa-tion between the hard and soft scales Λ( t ) (cid:29) ω pl ( t ) (withincreasing time t in our case) leads to an increasing sep-aration between plasmon mass and width, ω pl ( t ) (cid:29) γ ( t ).In the 2+1D cases, in contrast, the damping rate is ofthe order of the plasmon mass γ ( t ) ∼ ω pl ( t ) at all times.This indicates that the damping rate is not determinedby the hard degrees of freedom, but results from non-perturbative interactions between the soft gauge degreesof freedom. This is a qualitatively different dynamicalpicture than in three spatial dimensions and has the im-portant consequence that gluonic transverse and longitu-dinal excitations below the mass scale are too short-livedto form quasiparticles.Unlike the gauge fields, the scalar fields do have nar-row excitations at low momenta. At larger momenta,the damping rates saturate and the transverse and scalarexcitation peaks become identical.Similarly to gluonic plasmas in 3+1D, we observe thatthe statistical and spectral correlation functions obeya generalized fluctuation-dissipation relation. However,in contrast to 3+1D, the excitation peaks have a non-Lorentzian shape, that is reasonably well described byGaussian or hyperbolic secant distributions.We have also performed simulations in classical ther-mal equilibrium in 2+1D and observed qualitatively sim-ilar behaviour as along the self-similar attractor. We in-terpret this as a sign that the qualitative features maybe generic to 2+1D gauge theories.Our results indicate that an effective formulation of ki-netic theory in 2+1 dimensions should indeed be possiblefor large momenta p (cid:29) ω pl . However, this is complicatedby the fact that in lower dimensions infrared effects be-come more important. In particular, in 2+1D the dynam-ically generated screening scale gets equal contributionfrom all momenta, whereas in 3+1D the dominant con-tribution comes from hard particles. As a consequence,the HTL approximation breaks down already at the De-bye mass scale (instead of the magnetic scale in 3+1D).This picture is supported by our simulations. Therefore obtaining effective matrix elements for kinetic theory isa nonperturbative problem where the dynamics have tobe deduced in a self-consistent way taking into accountcontributions from the plasmon mass scale.We aim to return to the quasiparticle excitations of theGlasma in future work by measuring the spectral functionof expanding 2+1D systems. ACKNOWLEDGMENTS
We are grateful to A. Ipp, D. M¨uller, A. Rebhan andM. Strickland for valuable discussions. J. P. acknowl-edges the support by the International Office of the TUWien and would like to thank Institute for Theoreticalphysics for hospitality during part of this work. Thiswork has been supported by the European ResearchCouncil under grant no. ERC-2015-CoG-681707, by theEU Horizon 2020 research and innovation programme,STRONG-2020 project (grant agreement No 824093) andby the Academy of Finland, project 321840. This workwas funded in part by the Knut and Alice Wallenbergfoundation, contract number 2017.0036. The authorswish to acknowledge CSC - IT Center for Science, Fin-land, for computational resources. The computationalresults presented have been also achieved in part usingthe Vienna Scientific Cluster (VSC). The content of thisarticle does not reflect the official opinion of the Euro-pean Union and responsibility for the information andviews expressed therein lies entirely with the authors.
Appendix A: Formulas from the HTL framework1. Polarization tensor
In the HTL formalism, a crucial quantity is the po-larization tensor Π µν . For a non-Abelian SU( N ) gaugetheory in d spatial dimensions, its leading order (LO)expression reads [37, 38]Π µν ( K ) = g (cid:90) p V µ ∂ ¯ f p ∂ p β (cid:18) g νβ − V ν K β K ρ V ρ + i(cid:15) (cid:19) , (A1)where ∂/∂ p = 0, with Minkowski metric g νβ =diag(1 , − ), (cid:15) → V = (1 , p /ω p ), K = ( ω, k ) andwith the approximation ω p ≈ p for the dispersion re-lation, which is valid for sufficiently large momenta. Ingeneral, the frequency ω and momentum k are indepen-dent. The polarization tensor is gauge invariant, sym-metric Π µν = Π νµ and transverse with K µ Π µν ( K ) = 0for all ν . For the latter to be true, the distribution func-tion has to vanish at large momentum lim p i →∞ f p = 0for directions i , which is fulfilled in practice.We have summed different bosonic contributionsinto the distribution function ¯ f = N c (cid:80) λ f ( λ ) =: N c d pol f ( t, p ), where each f ( λ ) corresponds to the dis-tribution of the fields with non-longitudinal polarization4 λ . For the 2+1D theory, one has d pol = 1 and theintegration reads (cid:82) p = (cid:82) d p/ (2 π ) . When consideringthe Glasma-like theory, the distribution function has theform f ( t, p ) = 2 πδ ( p z ) f ( t, p ⊥ ), such that integrationbecomes identical to the 2+1D theory. Since the scalarcontributions correspond to a polarization in z direction,the distribution function f ( t, p ) becomes an average overgauge and scalar distributions, i.e., d pol f ( t, p ) = f G + f φ with d pol = 2 in the Glasma-like theory. Note that theresulting expressions are the same for both 2+1D andGlasma-like theories with only a different d pol .Using the isotropy of f ( t, p ) and performing an inte-gration by parts, the expression (A1) can be cast into theformΠ ij ( K ) = m D (cid:90) π d ϕ π (cid:20) δ ij − k i v j + k j v i − K ρ V ρ − i(cid:15) + ( − ω + k ) v i v j ( − K ρ V ρ − i(cid:15) ) (cid:21) , (A2)with i, j = 1 , m D = d pol N c (cid:90) d p (2 π ) g f ( t, p ) p . (A3)The transverse and longitudinal polarizations of Π ij ( K )can be computed asΠ T ( ω/k ) = q i Π ij ( K ) q j k , Π L ( ω/k ) = k i Π ij ( K ) k j ω , (A4)with the transverse vector q = ( − k , k ) such that q · k =0. Evaluating the residual integral [60], one arrives atΠ T ( x ) = m D x (cid:18) x − x − √ x + 1 √ x − (cid:19) (A5)Π L ( x ) = m D (cid:18) − x x − x + 1) / ( x − / (cid:19) , (A6)with x = ω/k .For the Glasma-like theory, we also need the scalarcomponent Π φ = Π zz . To obtain it, we start from (A1)and perform an integration by partsΠ φ = g (cid:90) d p (2 π ) ¯ f ( t, p ) (cid:20) ∂∂ p z v z + ∂∂ p l v z k l K ρ V ρ + i(cid:15) (cid:21) = m D = const , (A7)resulting from p z = k z = 0.
2. Retarded propagator
Next we consider the retarded propagator in temporalgauge G HTL ij ( ω, p ), as employed in this work, and ask howthese components of the polarization tensor emerge. Us-ing a tensor decomposition as suggested in Refs. [17, 61]with A ij = δ ij − p i p j p , B ij = p i p j p , C ij = n i n j n , (A8) where for the Glasma-like theory the normal vector issimply the unit vector in z direction n = e z , the propa-gator can be written as G HTL ( ω, p ) = G HTL T A + G HTL L B + G HTL φ C, (A9)with transverse, longitudinal and scalar polarizations G HTL T ( ω, p ) = − ω − p − Π T ( ω/p ) (A10) G HTL L ( ω, p ) = p ω − p − Π L ( ω/p ) (A11) G HTL φ ( ω, p ) = − ω − p − m D . (A12)One immediately observes the same screening propertiesin the static limit as in the 3 + 1-dimensional case, withlim ω → Π T = 0 , lim ω → − Π L = m D , (A13)additional to the static screening of scalar fields withlim ω → − ω + p + m D = p + m D . This gives the inter-pretation of m D as the Debye mass, as anticipated by itsnotation.
3. Spectral function
Each polarization of the spectral function can now becomputed as the imaginary part of the respective re-tarded propagator ρ HTL α ( ω, p ) = 2 Im G HTL α ( ω + i(cid:15), p ) , (A14)with α = T, L, φ denoting the polarization. They satisfythe sum rules˙ ρ T ( t, ∆ t =0 , p ) = (cid:90) ∞−∞ d ω π ωρ T ( t, ω, p ) = 1 (A15)˙ ρ L ( t, ∆ t =0 , p ) = (cid:90) ∞−∞ d ω π ωρ L ( t, ω, p ) = m D p + m D (A16)˙ ρ φ ( t, ∆ t =0 , p ) = (cid:90) ∞−∞ d ω π ωρ φ ( t, ω, p ) = 1 , (A17)which are the same as in 3+1D for T and L .The HTL spectral functions can be further split intoa Landau damping part for | ω | < p and a quasiparticlepart ρ HTL α ( ω, p ) = ρ Landau α ( ω, p )+ 2 πZ α ( p ) (cid:2) δ (cid:0) ω − ω HTL α ( p ) (cid:1) + δ (cid:0) ω + ω HTL α ( p ) (cid:1)(cid:3) , (A18) For the longitudinal spectral function, the correct (cid:15) prescriptionis to use Π L (( ω + i(cid:15) ) /p ) in G HTL L while keeping the prefactor p /ω in Eq. (A11) without this prescription. Otherwise onewould obtain a term proportional to δ ( ω ), which would lead towrong sum rules. ω HTL α ( p ) and residues Z α ( p )of quasiparticle excitations as discussed below.The Landau damping contributions have been writtenin the main text in Eqs. (20)-(22). The dispersion re-lations of quasiparticle excitations can be deduced from(A10)-(A12) by looking for poles in the propagators for ω > p . Defining ˜ p = p/m D , the dispersion relations canbe calculated analytically ω HTL T ( p ) = m D (cid:115) p − p + (cid:112) p − p (A19) ω HTL L ( p ) = m D p (cid:112) p (A20) ω HTL φ ( p ) = (cid:113) m D + p , (A21)where the dispersion relation for scalar particles is simplythe relativistic dispersion of free particles. It is interest-ing to study their behavior at low and high momenta.These are ω HTL T ( p ) p (cid:28) m D (cid:39) (cid:114) ω + 54 p (A22) ω HTL L ( p ) p (cid:28) m D (cid:39) (cid:114) ω + 34 p (A23) ω HTL T ( p ) p (cid:29) m D (cid:39) (cid:113) p + m (A24) ω HTL L ( p ) p (cid:29) m D (cid:39) (cid:115) p + m p , (A25)with ω = m D , m = m D . (A26)The leading asymptotic behavior is similar to the d = 3case. Like there, one has ω HTL T ( p ) > ω HTL L ( p ) for p > m D can also be interpreted as the asymptotic mass.On the other hand, the longitudinally polarized quasipar-ticles do not have an asymptotic mass, just as in the 3+1-dimensional case. An important difference to the latteris, however, that their dispersion relation approaches theultra-relativistic limit ω (cid:39) p considerably slower than inthe d = 3 case where the approach is exponential.The residues at the quasiparticle peaks are defined as Z α ( p ) = − (cid:20) ∂G − α ( ω, p ) ∂ω (cid:21) − ω = ω HTL α ( p ) (A27) and read Z T ( p ) = (cid:32) ω T − m D p (cid:32) ω T − ω T − p (cid:112) ω T − p (cid:33)(cid:33) − (A28) Z L ( p ) = (cid:0) ω L − p (cid:1) / ω L m D (A29) Z φ ( p ) = 12 ω φ , (A30)with asymptotic values Z T ( p ) p (cid:28) m D (cid:39) ω pl (A31) Z L ( p ) p (cid:28) m D (cid:39) ω pl (A32) Z T ( p ) p (cid:29) m D (cid:39) p (A33) Z L ( p ) p (cid:29) m D (cid:39) m D p , (A34)where we dropped the HTL label and the dependenceon momentum of the dispersion relations to shorten theexpressions. Thus, also in the 2+1D case one expectscontributions from longitudinal quasiparticles to decreaseat high momenta. However, it is a power law decrease,as opposed to the exponential decrease of 3+1D gaugetheory. Note that the leading contribution of Z T for p (cid:28) m D and p (cid:29) m D and Z L for p (cid:28) m D are thesame in 2+1D and 3+1D. This implies that the sum rules(A15)-(A16) for the spectral functions are dominated byquasiparticle contributions for p (cid:28) m D for both polar-izations and for p (cid:29) m D for transverse polarizations,while the longitudinal sum rule for p (cid:29) m D is domi-nated by the Landau damping contribution. ω / T C o rr e l a t i on ρ . T T Thermal, p = 0 T
Tt = 100Tt = 400 0 0.2 0.4 0.6 0.8 1
Thermal, p = 0.6 T
Tt = 100Tt = 400
FIG. 11. The correlator ˙ ρ T as a function of frequency for mo-menta p = 0 T (left) and p = 0 . T (right) for a 2+1D theoryin classical thermal equilibrium. Time translation invarianceis shown by comparing the correlations at two different times.One finds the same properties as for the non-equilibrium caseconsidered in the main part of the paper. Note that in theright panel the curve for T t = 100 is barely visible because itaccurately coincides with the curve for
T t = 400. Appendix B: Cross check: classical thermalequilibrium
Here we compute the correlators of the 2+1D the-ory in frequency space in classical thermal equilibrium,and compare the results to the non-equilibrium systemsdiscussed above. To obtain the classical thermal state,we initialize our system with f ( t = 0 , p ) = T /p on a256 lattice with lattice spacing T a = 0 . g = 0 . T . After the initialization we restore the Gausslaw constraint and let the system evolve. We find thatthe equal-time correlators (cid:104) EE (cid:105) T/L ( t, ∆ t =0 , p ) becomestationary for times T t (cid:38)
20, marking the onset of clas-sical thermal equilibrium. Due to the small Debye mass m D ∼ g T (cid:28) T , the resulting thermal state has almostthe same temperature as the initial temperature param-eter T ≈ T .Next, we extract correlation functions at differenttimes T t = 100 and 400 to demonstrate time-translationinvariance explicitly, as should be the case in thermalequilibrium. Our results are shown for the transverse dot-ted spectral function ˙ ρ T ( t, ω, p ) in Fig. 11 for momenta p = 0 (left) and p = 0 . T (right). One observes thatthe curves corresponding to different extraction times lieindeed on top of each other within statistical uncertain-ties. For higher momenta, where error bars are small, thecurves can be barely distinguished by eye (right panel).This confirms that we have reached classical thermal equilibrium.Our numerical results qualitatively agree with our find-ings in the non-equilibrium systems above. We can sum-marize our most important findings as follows:1. We observe a fluctuation-dissipation relation, asshould be trivially valid in thermal equilibrium.2. The gluonic excitations are broad and their widthis of the same order as the plasmon mass.3. The dotted spectral function is seen to approach afinite value for ω → [1] F. Gelis, E. Iancu, J. Jalilian-Marian andR. Venugopalan, The color glass condensate , Ann. Rev.Nucl. Part. Sci. (2010) 463 [ arXiv:1002.0333[hep-ph] ].[2] T. Lappi and L. McLerran, Some features of the glasma , Nucl. Phys.
A772 (2006) 200 [ arXiv:hep-ph/0602189 ].[3] F. Gelis,
Initial state and thermalization in the colorglass condensate framework , Int. J. Mod. Phys.
E24 (2015) 1530008 [ arXiv:1508.07974 [hep-ph] ].[4] A. Krasnitz and R. Venugopalan,
Non-perturbativecomputation of gluon mini-jet production in nuclearcollisions at very high energies , Nucl. Phys.
B557 (1999) 237 [ arXiv:hep-ph/9809433 ].[5] A. Krasnitz, Y. Nara and R. Venugopalan,
Coherentgluon production in very high energy heavy ioncollisions , Phys. Rev. Lett. (2001) 192302[ arXiv:hep-ph/0108092 ].[6] T. Lappi, Production of gluons in the classical fieldmodel for heavy ion collisions , Phys. Rev.
C67 (2003)054903 [ arXiv:hep-ph/0303076 ].[7] A. Krasnitz, Y. Nara and R. Venugopalan,
Classicalgluodynamics of high energy nuclear collisions: Anerratum and an update , Nucl. Phys.
A727 (2003) 427[ arXiv:hep-ph/0305112 ].[8] T. Lappi,
Energy density of the glasma , Phys. Lett.
B643 (2006) 11 [ arXiv:hep-ph/0606207 ].[9] R. J. Fries, J. I. Kapusta and Y. Li,
Near-fields andinitial energy density in the color glass condensatemodel , arXiv:nucl-th/0604054 . [10] G. Chen, R. J. Fries, J. I. Kapusta and Y. Li, Earlytime dynamics of gluon fields in high energy nuclearcollisions , Phys. Rev.
C92 (2015) 064912[ arXiv:1507.03524 [nucl-th] ].[11] B. Schenke, P. Tribedy and R. Venugopalan,
Fluctuating glasma initial conditions and flow in heavyion collisions , Phys. Rev. Lett. (2012) 252301[ arXiv:1202.6646 [nucl-th] ].[12] B. Schenke, P. Tribedy and R. Venugopalan,
Event-by-event gluon multiplicity, energy density, andeccentricities in ultrarelativistic heavy-ion collisions , Phys. Rev.
C86 (2012) 034908 [ arXiv:1206.6805[hep-ph] ].[13] H. M¨antysaari, B. Schenke, C. Shen and P. Tribedy,
Imprints of fluctuating proton shapes on flow inproton-lead collisions at the LHC , Phys. Lett.
B772 (2017) 681 [ arXiv:1705.03177 [nucl-th] ].[14] D. Gelfand, A. Ipp and D. M¨uller,
Simulating collisionsof thick nuclei in the color glass condensate framework , Phys. Rev.
D94 (2016) 014020 [ arXiv:1605.07184[hep-ph] ].[15] A. Ipp and D. M¨uller,
Broken boost invariance in theglasma via finite nuclei thickness , Phys. Lett.
B771 (2017) 74 [ arXiv:1703.00017 [hep-ph] ].[16] A. Ipp and D. M¨uller,
Implicit schemes for real-timelattice gauge theory , Eur. Phys. J.
C78 (2018) 884[ arXiv:1804.01995 [hep-lat] ].[17] P. Romatschke and M. Strickland,
Collective modes ofan anisotropic quark-gluon plasma , Phys. Rev.
D68 (2003) 036004 [ arXiv:hep-ph/0304092 ].[18] A. Rebhan, P. Romatschke and M. Strickland, Hard-loop dynamics of non-Abelian plasma instabilities , Phys. Rev. Lett. (2005) 102303[ arXiv:hep-ph/0412016 [hep-ph] ].[19] P. Romatschke and R. Venugopalan, Collectivenon-Abelian instabilities in a melting color glasscondensate , Phys. Rev. Lett. (2006) 062302[ arXiv:hep-ph/0510121 ].[20] T. Epelbaum and F. Gelis, Fluctuations of the initialcolor fields in high energy heavy ion collisions , Phys.Rev.
D88 (2013) 085015 [ arXiv:1307.1765 [hep-ph] ].[21] T. Epelbaum and F. Gelis,
Pressure isotropization inhigh energy heavy ion collisions , Phys. Rev. Lett. (2013) 232301 [ arXiv:1307.2214 [hep-ph] ].[22] A. Mueller and D. Son,
On the equivalence between theBoltzmann equation and classical field theory at largeoccupation numbers , Phys. Lett. B (2004) 279[ arXiv:hep-ph/0212198 ].[23] S. Jeon,
The Boltzmann equation in classical andquantum field theory , Phys. Rev.
C72 (2005) 014907[ arXiv:hep-ph/0412121 [hep-ph] ].[24] J. Berges,
Introduction to nonequilibrium quantum fieldtheory , AIP Conf. Proc. (2005) 3[ arXiv:hep-ph/0409233 [hep-ph] ].[25] A. Kurkela and G. D. Moore,
UV cascade in classicalYang-Mills theory , Phys. Rev.
D86 (2012) 056008[ arXiv:1207.1663 [hep-ph] ].[26] J. Berges, S. Scheffler and D. Sexty,
Bottom-upisotropization in classical-statistical lattice gauge theory , Phys. Rev.
D77 (2008) 034504 [ arXiv:0712.3514[hep-ph] ].[27] J. Berges, S. Schlichting and D. Sexty,
Over-populatedgauge fields on the lattice , Phys. Rev.
D86 (2012)074006 [ arXiv:1203.4646 [hep-ph] ].[28] S. Schlichting,
Turbulent thermalization of weaklycoupled non-abelian plasmas , Phys. Rev.
D86 (2012)065008 [ arXiv:1207.1450 [hep-ph] ].[29] J. Berges, K. Boguslavski, S. Schlichting andR. Venugopalan,
Universal attractor in a highlyoccupied non Abelian plasma , Phys. Rev.
D89 (2014)114007 [ arXiv:1311.3005 [hep-ph] ].[30] J. Berges, K. Boguslavski, S. Schlichting andR. Venugopalan,
Turbulent thermalization process inheavy-ion collisions at ultrarelativistic energies , Phys.Rev.
D89 (2014) 074011 [ arXiv:1303.5650 [hep-ph] ].[31] J. Berges, K. Boguslavski, S. Schlichting andR. Venugopalan,
Basin of attraction for turbulentthermalization and the range of validity ofclassical-statistical simulations , JHEP (2014) 054[ arXiv:1312.5216 [hep-ph] ].[32] A. Kurkela and Y. Zhu, Isotropization andhydrodynamization in weakly coupled heavy-ioncollisions , Phys. Rev. Lett. (2015) 182301[ arXiv:1506.06647 [hep-ph] ].[33] A. Kurkela, A. Mazeliauskas, J.-F. Paquet,S. Schlichting and D. Teaney,
Matching thenonequilibrium initial stage of heavy ion collisions tohydrodynamics with QCD kinetic theory , Phys. Rev.Lett. (2019) 122302 [ arXiv:1805.01604 [hep-ph] ].[34] A. Kurkela, A. Mazeliauskas, J.-F. Paquet,S. Schlichting and D. Teaney,
Effective kineticdescription of event-by-event pre-equilibrium dynamicsin high-energy heavy-ion collisions , Phys. Rev.
C99 (2019) 034910 [ arXiv:1805.00961 [hep-ph] ].[35] R. Baier, A. H. Mueller, D. Schiff and D. Son, ‘Bottomup’ thermalization in heavy ion collisions , Phys. Lett. B (2001) 51 [ arXiv:hep-ph/0009237 ].[36] K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron,
Highly occupied gauge theories in 2+1 dimensions: Aself-similar attractor , Phys. Rev.
D100 (2019) 094022[ arXiv:1907.05892 [hep-ph] ].[37] E. Braaten and R. D. Pisarski,
Soft amplitudes in hotgauge theories: A general analysis , Nucl. Phys.
B337 (1990) 569.[38] J.-P. Blaizot and E. Iancu,
The quark gluon plasma:Collective dynamics and hard thermal loops , Phys. Rept. (2002) 355 [ arXiv:hep-ph/0101103 [hep-ph] ].[39] A. Kurkela, T. Lappi and J. Peuron,
Time evolution oflinearized gauge field fluctuations on a real-time lattice , Eur. Phys. J.
C76 (2016) 688 [ arXiv:1610.01355[hep-lat] ].[40] K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron,
Spectral function for overoccupied gluodynamics fromreal-time lattice simulations , Phys. Rev.
D98 (2018)014006 [ arXiv:1804.01966 [hep-ph] ].[41] A. Pi˜neiro Orioli and J. Berges,
Breaking thefluctuation-dissipation relation by universal transportprocesses , Phys. Rev. Lett. (2019) 150401[ arXiv:1810.12392 [cond-mat.quant-gas] ].[42] K. Boguslavski and A. Pi˜neiro Orioli,
Unraveling thenature of universal dynamics in O ( N ) theories , Phys.Rev. D (2020) 091902 [ arXiv:1911.04506[hep-ph] ].[43] G. Aarts,
Spectral function at high temperature in theclassical approximation , Phys. Lett.
B518 (2001) 315[ arXiv:hep-ph/0108125 [hep-ph] ].[44] S. Schlichting, D. Smith and L. von Smekal,
Spectralfunctions and critical dynamics of the O (4) model fromclassical-statistical lattice simulations , Nucl. Phys. B (2020) 114868 [ arXiv:1908.00912 [hep-lat] ].[45] K. Boguslavski, A. Kurkela, T. Lappi and J. Peuron,
Heavy quark diffusion in an overoccupied gluon plasma , JHEP (2020) 077 [ arXiv:2005.02418 [hep-ph] ].[46] J. Berges, A. Rothkopf and J. Schmidt, Non-thermalfixed points: Effective weak-coupling for stronglycorrelated systems far from equilibrium , Phys. Rev. Lett. (2008) 041603 [ arXiv:0803.0131 [hep-ph] ].[47] A. Pi˜neiro Orioli, K. Boguslavski and J. Berges,
Universal self-similar dynamics of relativistic andnonrelativistic field theories near nonthermal fixedpoints , Phys. Rev. D (2015) 025041[ arXiv:1503.02498 [hep-ph] ].[48] R. Walz, K. Boguslavski and J. Berges, Large-N kinetictheory for highly occupied systems , Phys. Rev. D (2018) 116011 [ arXiv:1710.11146 [hep-ph] ].[49] J. Berges, K. Boguslavski, A. Chatrchyan andJ. Jaeckel, Attractive versus repulsive interactions in theBose-Einstein condensation dynamics of relativistic fieldtheories , Phys. Rev. D (2017) 076020[ arXiv:1707.07696 [hep-ph] ].[50] I. Chantesana, A. Pi˜neiro Orioli and T. Gasenzer, Kinetic theory of nonthermal fixed points in a Bose gas , Phys. Rev. A (2019) 043620 [ arXiv:1801.09490[cond-mat.quant-gas] ].[51] J. Berges, K. Boguslavski, S. Schlichting andR. Venugopalan, Universality far from equilibrium:From superfluid Bose gases to heavy-ion collisions , Phys. Rev. Lett. (2015) 061601 [ arXiv:1408.1670[hep-ph] ].[52] J. Berges, K. Boguslavski, S. Schlichting andR. Venugopalan,
Nonequilibrium fixed points inlongitudinally expanding scalar theories: infraredcascade, Bose condensation and a challenge for kinetictheory , Phys. Rev.
D92 (2015) 096006[ arXiv:1508.03073 [hep-ph] ].[53] M. Pr¨ufer, P. Kunkel, H. Strobel, S. Lannig,D. Linnemann, C.-M. Schmied, J. Berges, T. Gasenzerand M. K. Oberthaler,
Observation of universaldynamics in a spinor Bose gas far from equilibrium , Nature (2018) 217 [ arXiv:1805.11881[cond-mat.quant-gas] ].[54] S. Erne, R. B¨ucker, T. Gasenzer, J. Berges andJ. Schmiedmayer,
Universal dynamics in an isolatedone-dimensional Bose gas far from equilibrium , Nature (2018) 225 [ arXiv:1805.12310[cond-mat.quant-gas] ].[55] J. A. Glidden, C. Eigen, L. H. Dogra, T. A. Hilker,R. P. Smith and Z. Hadzibabic,
Bidirectional dynamicscaling in an isolated Bose gas far from equilibrium , arXiv:2006.01118 [cond-mat.quant-gas] .[56] T. Lappi and J. Peuron, Plasmon mass scale in classicalnonequilibrium gauge theory , Phys. Rev.
D95 (2017)014025 [ arXiv:1610.03711 [hep-ph] ].[57] T. Lappi and J. Peuron,
Plasmon mass scale in twodimensional classical nonequilibrium gauge theory , Phys. Rev. D (2018) 034017 [ arXiv:1712.02194[hep-lat] ].[58] L. Shen and J. Berges, Spectral, statistical and vertexfunctions in scalar quantum field theory far fromequilibrium , Phys. Rev. D (2020) 056009[ arXiv:1912.07565 [hep-ph] ].[59] P. B. Arnold, G. D. Moore and L. G. Yaffe,
Effectivekinetic theory for high temperature gauge theories , JHEP (2003) 030 [ arXiv:hep-ph/0209353 [hep-ph] ].[60] P. Romatschke and M. Strickland, Collective modes ofan anisotropic quark-gluon plasma. II , Phys. Rev.
D70 (2004) 116006 [ arXiv:hep-ph/0406188 ].[61] R. Kobes, G. Kunstatter and A. Rebhan,
Gaugedependence identities and their application at finitetemperature , Nucl. Phys.