Hydrodynamization and non-equilibrium Green's functions in kinetic theory
Syo Kamata, Mauricio Martinez, Philip Plaschke, Stephan Ochsenfeld, S. Schlichting
HHydrodynamization and non-equilibrium Green’s functions in kinetic theory
Syo Kamata,
1, 2, ∗ Mauricio Martinez, † Philip Plaschke, ‡ Stephan Ochsenfeld, § and S¨oren Schlichting ¶ College of Physics and Communication Electronics,Jiangxi Normal University, Nanchang 330022, China Department of Physics, North Carolina State University, Raleigh, NC 27695, USA Fakult¨at f¨ur Physik, Universit¨at Bielefeld, D-33615 Bielefeld, Germany (Dated: April 16, 2020)Non-equilibrium Green’s functions provide an efficient way to describe the evolution of the energy-momentum tensor during the early time pre-equilibrium stage of high-energy heavy ion collisions.Besides their practical relevance they also provide a meaningful way to address the question when andto what extent a hydrodynamic description of the system becomes applicable. Within the kinetictheory framework we derive a new method to calculate time dependent non-equilibrium Green’sfunctions describing the evolution of energy and momentum perturbations on top of an evolvingfar-from-equilibrium background. We discuss the approach towards viscous hydrodynamics alongwith the emergence of various scaling phenomena for conformal systems. By comparing our resultsobtained in the relaxation time approximation to previous calculations in Yang-Mills kinetic theory,we further address the question which macroscopic features of the energy momentum tensor aresensitive to the underlying microscopic dynamics.
I. INTRODUCTION
Ultrarelativistic heavy ion experiments carried out atthe Relativistic Heavy Ion Collider (RHIC) and the LargeHadron Collider (LHC) have produced a new state ofmatter where quarks and gluons are liberated from theincoming nuclei [1–4] Since the lifetime of this Quark-Gluon Plasma (QGP) is very short, its properties arereconstructed by analyzing a large set of hadronic ob-servables. The phenomenological studies of the collectedwealth of experimental data have shown that hydrody-namical models provide robust tools to explain the dy-namics of the QGP [5–9]. As a result, a new paradigmin high energy nuclear physics has emerged where hydro-dynamics plays a central role.One of the most consistent findings in the hydrody-namical modeling of heavy ion collisions is a small valueof the shear viscosity over entropy ratio η/s ∼ / (4 π ),extracted from large number of observables measured inultrarelativistic heavy ion collisions over a different rangeof energies at both RHIC and LHC. Currently, the stan-dard approach to describe high-energy heavy ion colli-sions is based on multistage models where the space-timedynamics of the QGP is described using relativistic vis-cous hydrodynamics. Subsequently, the energy density ofthe QGP is converted into a hadron gas which eventuallyfreezes out, yielding the final state hadronic observablesthat can compared to experimental measurements. Onekey ingredient in this phenomenological approach are theinitial conditions, characterizing the initial distributionsof the energy density and flow velocities, which are then ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] propagated via hydrodynamics. Clearly, the calculationof these quantities from first principles QCD representsan enormous challenge, which despite important develop-ments during the last years [10–18] has not been answeredcompletely.One important feature shared by all initial state mod-els is related to the far-from-equilibrium nature of theQCD matter, produced immediately after the collisionof heavy nuclei on a time scale τ (cid:28) / c [19]. Sinceat these early times the QCD matter is also subject toa rapid longitudinal expansion, viscous hydrodynamicswhich is an effective theory for the long-time and longwave-length behavior close to equilibrium is not neces-sarily applicable at such early times. Although in recentyears the existence of a far-from-equilibrium fluid dynam-ical theory has been advocated [20–42], its formulationremains to be completed and thus, in practice hydrody-namic simulations of heavy-ion collisions are usually ini-tialized after a certain period time τ = τ hydro ∼ / c,where the long. expansion is less rapid and the QGP hasevolved towards local thermal equilibrium.The question how to match the non-equilibratium ini-tial state at τ to the initial state for hydrodynamics at τ hydro was the subject of a series of recent papers [43–46]. Based on an underlying microscopic descriptionin QCD kinetic theory, the initial conditions for theenergy-momentum tensor T µν ( τ hydro ) at the beginningof the hydrodynamic phase are hereby obtained fromthe initial energy-momentum tensor T µν ( τ ), by deter-mining the evolution of macroscopic quantities, based on non-equilibrium Green’s functions of the energy momen-tum tensor [44, 45]. This new framework, dubbed asKøMPøST [44, 45], has proven to be a powerful toolto describe the pre-equilibrium evolution of heavy-ioncollisions on an event-by-event basis [47, 48]. Whilethe detailed phenomenological consequences of the pre-equilibrium stage remain to be explored, it has beendemonstrated that the subsequent hydrodynamic evolu- a r X i v : . [ h e p - ph ] A p r tion becomes independent of the hydrodynamic initial-ization time τ hydro as long as the latter is chosen to be inthe regime where both kinetic theory and hydrodynamicdescriptions overlap.While [44, 45] presented a calculation of non-equilibrium Green’s functions in QCD kinetic theory, it isalso interesting to explore to what extent the microscopicdetails affect the evolution of the macroscopic quantitiesfar-from-equilibrium. In this work we therefore present anew method to calculate non-equilibrium Green’s func-tions of the energy momentum tensor. We generalizethe method of moments approach [49, 50] by incorpo-rating the response of a far-from-equilibrium expandingplasma against linear perturbations, and analyze theirnon-equilibrium in the relaxation time approximation.Instead of solving linear kinetic theory for perturbationsof the phase-space distribution function we analyze theequations of motion of the corresponding linearized mo-ments. Eventually, the Green’s functions of the energymomentum tensor are reconstructed from the linearizedmoments.This paper is organized as follows: In Sect. II we in-troduce the relevant aspects of the Boltzmann equation.We explain the main aspects of the emergent attract-ing behavior of the non-equilibrated Bjorken flow back-ground in Sect. III. In Sect. IV we present the formalismto describe the space-time evolution of the perturbations.Using this formalism we proceed in Sect. V to calculatethe Green’s functions of the energy-momentum tensor.Conclusions and outlook are discussed in Sect. VI. Someof technical aspects of our work are briefly discussed inAppendix A. II. BOLTZMANN EQUATION WITHIN THERELAXATION TIME APPROXIMATION
Starting point of our analysis is the Boltzmann equa-tion within the relaxation time approximation (RTA) p µ ∂ µ f = C [ f ] = − p µ u µ ( x ) τ R (cid:104) f − f eq ( p µ β µ ( x )) (cid:105) , (1)where the coordinate system defined in Minkowskispace is x µ = ( x , x , x ) with the metric g µν =diag(1 , − , − , − β µ ( x ) = u µ ( x ) /T ( x ) with the local rest-frame velocity u µ ( x ) de-termined via the Landau matching condition T µν ( x ) u ν ( x ) = e ( x ) u µ ( x ) , (2a) e ( x ) = e eq ( T ( x )) . (2b)The fluid velocity is defined as a time-like eigenvector( u = +1) of the energy momentum tensor for on-shellparticles T µν ( x ) = (cid:104) p µ p ν (cid:105) , where we denote the on-shell momentum average of anyobservable as (cid:104) O (cid:105) X = ν eff (2 π ) (cid:90) d p (cid:112) − g ( x ) (2 π ) δ ( p ) 2 θ ( p ) × O ( x µ , p µ ) f X ( x µ , p µ ) . (3)The effective temperature T ( x ) entering in Eq. (2b) is de-termined from the (equilibrium) equation of state e eq ( T ),by matching the corresponding eigenvalue e to the equi-librium energy density. If not stated otherwise, we willconsider an ultra-relativistic system of massless bosonswhere the equilibrium distribution function is given bythe Bose-Einstein distribution f eq. ( x ) = 1 / ( e x − e eq ( T ) = ν eff π T . (4)We are interested in longitudinal boost-invariant expand-ing system and thus, we use the hyperbolic Bjorken co-ordinates defined in terms of the cartesian coordinatesas τ = (cid:113) ( x ) − ( x ) ς = arctanh (cid:0) x /x (cid:1) , (5)so the metric g µν = diag (1 , − , − , − τ ) and (cid:112) − g ( x ) = τ . Similarly, the four-momentum of a relativistic mass-less particle is p µ = ( p T cosh( y ) , p , p T sinh( y )) (6)with y = arctanh ( p /p ) and p T = | p | .In the Bjorken coordinates the Boltzmann equation (1)takes the following form (cid:2) p τ ∂ τ + p i ∂ i + p ς ∂ ς (cid:3) f ( x, p ) = (7) − p µ u µ ( x ) τ R (cid:104) f ( x, p ) − f eq ( p µ β µ ( x )) (cid:105) . where p τ = p T cosh( y − ς ) , p ς = 1 τ p T sinh( y − ς ) (8)Hereafter we shall use the roman letter i = x, y to denotethe summation only over the transverse coordinates. Byvirtue of the coordinate transformation (5), the deriva-tives w.r.t. τ, ς in Eq. (7) are taken at constant p T andmomentum space rapidity y . However, when analyzingthe dynamics of a boost invariant medium it is moreconvenient to work with the (dimensionless) longitudi-nal momentum variable p ς = − τ p ς = − τ p T sinh( y − ς ) . (9)By transforming the space-time derivatives according to ∂ τ = ∂ τ | p ς − p T sinh( y − ς ) ∂ p ς , (10) ∂ ς = ∂ ς | p ς + τ p T cosh( y − ς ) ∂ p ς , (11)one finds that the two additional terms cancel each other,such that (cid:104) p τ ∂ τ + p i ∂ i − p ς τ ∂ ς (cid:105) f ( x, p ) = (12) − p µ u µ ( x ) τ R (cid:104) f ( x, p ) − f eq ( p µ β µ ( x )) (cid:105) , which is the form of the equation that we will considerin this work. We note in passing, that it is also commonin the literature to express the dynamics in terms of thelongitudinal momentum in the local rest frame p (cid:107) = τ p ς = p T sinh( y − ς ) . (13)In this case only the derivative w.r.t. the longitudinalrapidity is affected by the transformation ∂ ς = ∂ ς | p (cid:107) − p T cosh( y − ς ) ∂ p (cid:107) , (14)and the Boltzmann equation takes the form (see e.g. [51]) (cid:20) p τ ∂ τ + p i ∂ i + p (cid:107) τ ∂ ς − p τ p (cid:107) τ ∂ p (cid:107) (cid:21) f ( x, p ) = (15) − p µ u µ ( x ) τ R (cid:104) f ( x, p ) − f eq ( p µ β µ ( x )) (cid:105) . III. EVOLUTION OF BOOST INVARIANTHOMOGENOUS BACKGROUND
During the early pre-equilibrium stage of a heavy-ioncollision, which lasts about 1fm / c, the non-equilibriumplasma is subject to a rapid longitudinal expansion. Con-versely, in the transverse plane the plasma is initiallycreated at rest, and the transverse expansion only buildsup in response to local energy-density gradients on a timescale τ ∼ R , where R denotes the system size. Due to thisseparation of scales, it is a reasonable assumption to ne-glect the transverse expansion during the pre-equilibriumstage, and first consider an idealized situation of Bjorkenflow, where the system is longitudinally boost-invariant,parity invariant under spatial reflexions along the longi-tudinal beam line and azimuthally symmetric and trans-lationally invariant in the transverse plane. Space-timedependent variations can subsequently be addressed bystudying small deviations from this average backgroundbehavior, and will be discussed in Sec. IV.In this section we discuss the space-time evolution ofan expanding background undergoing Bjorken expansion,such that the aforementioned symmetries constrain thefunctional form of the distribution function for the back-ground to be of the form f ( x, p ) = f BG ( τ, p T , | p ς | ) . (16)The energy momentum tensor has only non-vanishing di-agonal components, i.e., T µνBG = diag ( e, p T , p T , p L /τ ) , (17) with the energy density ( e ), transverse and longitudinalpressures ( p T/L ) determined by e = T ττBG = (cid:10) ( p τ ) (cid:11) BG , (18a) p T = T xxBG = T yyBG = 12 (cid:10) ( p T ) (cid:11) BG , (18b) p L = T ςςBG = (cid:28) (cid:16) p ς τ (cid:17) (cid:29) BG . (18c)Due to scale invariance, the previous expressions auto-matically satisfy the tracelessness condition e = 2 p T + p L .Based on the explicit form of T µνBG the Landau matchingcondition in Eq. (2) becomes trivial with u µ = ( u τ , u i , u ς ) = (1 , , , , (19a) e = e eq , (19b)and the kinetic equation for the evolution of the back-ground distribution takes the familiar form τ ∂ τ f BG ( τ, p T , | p ς | ) = (20) − ττ R (cid:20) f BG (cid:16) τ, p T , | p ς | (cid:17) − f eq (cid:16) p τ T ( τ ) (cid:17)(cid:21) . where p τ = (cid:112) p T + ( p ς /τ ) according to the on-shellmass condition. A. Evolution equations for moments
Several strategies have been explored in the literatureto solve Eq. (20) [52]. Here we follow Grad’s originalapproach [50] where instead of finding solutions for thephase-space distribution function f ( x, p ) we study thedynamics of its moments. The latter are directly con-nected to macroscopic observables as we indicate below.We consider the following moments C ml ( τ ) = ν eff (cid:90) dp ς (2 π ) (cid:90) d p (2 π ) τ / (cid:113) p T + ( p ς /τ ) Y ml ( φ p , θ p ) f BG ( τ, p T , | p ς | ) , (21)where the angles are defined as cos θ p = p ς / ( τ p τ ),tan φ p = p /p in the co-moving coordinates and Y ml denote the spherical harmonics Y ml ( φ, θ ) = y ml P ml (cos( θ )) e imφ , (22)with normalization y ml = (cid:115) (2 l + 1)( l − m )!4 π ( l + m )! . (23)Non-vanishing components of the background energy-momentum tensor in (17) are related with the moments C ml in (21) as follows C ( τ ) = (cid:114) π τ / e ( τ ) , (24) C ( τ ) = (cid:114) π τ / (cid:104) p L ( τ ) − e ( τ ) (cid:105) . (25)Specifcially, if the background distribution function is inthermal equilibrium, one has C ml | eq ( τ ) = τ / e ( τ ) √ π δ l δ m , (26)where the Landau matching conditions (2) was enforced.Now the evolution equations of those moments is simplyobtained by taking the explicit time derivative in its def-inition. After some careful algebra and using a series ofidentities of the spherical harmonics (see App. A), theevolution equation for the moments C ml takes the follow-ing form τ ∂ τ C ml = b ml, − C ml − + b ml, C ml + b ml, +2 C ml +2 − ττ R (cid:16) C ml − C ml | eq (cid:17) , (27)where the coefficients b ml, − , b ml, and b ml, +2 are given by b ml, − = l + 22 l + 1 (cid:115) Γ( l − m + 1)Γ( l + m + 1)(2 l + 1)(2 l − l − m − l + m − b ml, = − l ( l + 1) − m l ( l + 1) − b ml, +2 = − l − l + 3 (cid:115) Γ( l − m + 3)Γ( l + m + 3)(2 l + 1)(2 l + 5)Γ( l − m + 1)Γ( l + m + 1)We note that in Eq. (27) only the moments l and l ± m arenot coupled and only moments with m = 0 are non-vanishing. Furthermore, Eqs. (27) present interestingmathematical resurgent properties which were discussedextensively in Refs. [53–55]. B. Initial conditions
Since at early times τ (cid:28) τ R the system is unable tosustain sizeable longitudinal momenta, the physically rel-evant initial conditions for the phase-space distributionare naturally of a form where p (cid:107) (cid:28) p T . Indeed, previ-ous works [35–37, 53, 55–57] have shown that the extremelimit where the (longitudinal) support of the phase-spacedistribution shrinks to a Dirac delta function correspondsto a non-equilibrium attractor of the kinetic equation.We will therefore consider precisely this case and choosethe initial phase-space distribution as f BG ( τ , p T , p ς ) = (2 π ) ν eff δ ( p ς ) dN dςd p d x , (29)where the normalization chosen such that the initial en-ergy density per unit rapidity remains constant, i.e. dE dςd x = lim τ → τ e ( τ ) = ( τ e ) = const . (30) By construction, the initial condition in Eq. (29) fixesthe initial values of the moments C ml at τ , i.e., C ml ( τ ) = τ / ( τ e ) y ml P ml (0) δ m (31)where P ml (0) = 2 m √ π Γ (cid:0) − l − m (cid:1) Γ (cid:0) − m − l (cid:1) , (32)such that for m = 0 one has P l (0) = sin( π/ l + 1)) √ π Γ(1 / l/ l/ . (33)Interestingly, one finds that at the level of relaxation timeapproximation, the shape of the azimuthally symmetricmomentum distribution dNdςd p d x is completely irrelevantfor the dynamics. Noteably, this is in sharp contrast tomore realistic description in Yang-Mills kinetic theory[23, 44, 45, 58], where different processes (e.g. inelasticand elastic interactions) exhibit different parametricdependencies on the momenta (see e.g. [19, 59]). C. Evolution of energy momentum tensor forconstant relaxation time ( τ R = const ) We first analyze the evolution of the energy momen-tum tensor for a constant relaxation time τ R = const . Bytruncating the evolution equations up to a certain valueof l < l max , setting C ml = C ml | eq = 0, we can obtain a nu-merical solution of the coupled set of evolution equationas a function of the dimensionless evolution time variable τ /τ R . Noteably the truncation scheme converges rapidlyexcept for very small values of τ /τ R (cid:28) l ≤ l max (cid:38)
32, and ifnot stated otherwise we use l max = 512 to produce thefigures.
1. Comparison with hydrodynamics
Our results are compactly summarized in Fig. 1, whichshows the non-equilibrium attractor solution for the evo-lution of the energy and pressure densities. For theBjorken flow we can also compare the numerical solutionof the Boltzmann equation in relaxation time approxi-mation to the truncation in relativistic viscous hydrody-namics. Starting from the conservation equation τ ∂ τ e = − e − p L , (34)the longitudinal and transverse pressure within the hy-drodynamic approach are given by p L = p + π ςς , p T = p − π ςς T e3p L E ne r g y / p r e ss u r e den s i t y : τ / T µ ν / ( e τ / ) ∞ Evolution time: x= τ / τ R Boltzmann RTA ( τ R =const)Free Streaming -- C -1 ∞ x Navier-Stokes Hydro -- 1-16/45x 0 0.2 0.4 0.6 0.8 1 1.2 0.01 0.1 1 10 E ff. r e l a x a t i on r a t e : τ R e ff / τ R = / τ / τ R ( p T - p L ) / ( e + p ) Evolution time: x= τ / τ R Boltzmann RTA ( τ R =const)l max =4l max =8Navier-Stokes Hydro FIG. 1. (top) Evolution of energy-momentum tensor forconst. relaxation time. (bottom) Effective relaxation timedetermined from the appropriate ratio of inverse Reynoldsand the Knudsen number (see text). where π ςς is the only independent component of the shearviscous tensor. Now π ςς is expanded up to second orderin the gradient hydrodynamical expansion as [60] π ςς = − ητ + 89 1 τ ( λ − ητ π ) , (36)where for the system under consideration the equation ofstate and transport coefficients in the RTA approxima-tion are determined by [37, 54, 61, 62] p = e/ , η/s = τ R T , T s = e + p ,τ π = τ R , λ = 335 τ R T s (37)such that η/τ = τ R τ e . Hence the asymptotic solu-tion for the energy density up to second order within thegradient hydrodynamic expansion (or alternatively when τ /τ R (cid:29)
1) takes the form τ / e ( τ ) ≈ (cid:16) τ / e (cid:17) ∞ (cid:18) − τ R τ − (cid:16) τ R τ (cid:17) (cid:19) , (38)with the asymptotic integration constant (cid:0) τ / e (cid:1) ∞ de-termined aslim τ →∞ τ / e ( τ ) = (cid:16) τ / e (cid:17) ∞ (cid:39) . τ / R ( τ e ) (39)from the numerical solution of the Boltzmann equationin relaxation time approximation. By comparing the dif-ferent curves in Fig. 1 one observes that viscous hydro-dynamics starts to describe the evolution of the energymomentum tensor around evolution times τ /τ R ∼ − τ R /τ in Eq. (38).Noteably, one can also define an effective ratio betweenthe inverse Reynolds and the Knudsen number as follows Re − Kn = 4 ττ R (cid:18) | π ςς | e + p (cid:19) , (40)with Kn = τ R /τ being the Knudsen number and theinverse Reynolds number Re − is Re − ∼ | π ςς | p = 23 | p L − p T | p . (41)Near equilibrium where the gradient expansion (36)holds, one can express the inverse Reynolds number as Re − = 163 ηs Kn + O ( Kn ) . (42)While the gradient expansion (36) breaks down at earlytimes, where the Knudsen number Kn (cid:29) τ effR away from equilibrium τ effR τ R = 1516 (cid:12)(cid:12)(cid:12)(cid:12) Re − Kn (cid:12)(cid:12)(cid:12)(cid:12) . (43)which is constructed such that at late times the ratio τ effR τ R in Eq. (43) converges to unity as expected.Numerical results for τ effR τ R shown in the bottom panelof Fig. 1 indicate that at early times where the Knudsennumber Kn (cid:29)
1, the effective relaxation rate is signifi-cantly reduced by the inclusion of higher order dynamicalmoments C ml . By explicitly comparing different trunca-tions ( l max = 4 , , · · · ) of the infinite hierarchy of mo-ment equations, one also observes a rapid convergence inthe sense that the evolution of the low order momentsbecomes increasingly insensitive to the higher order mo-ments. T e3p L E ne r g y / p r e ss u r e den s i t y : τ / T µ ν / ( e τ / ) ∞ Evolution time: w~=T( τ ) τ /(4 π η /s)Boltzmann RTA (T τ R =const)KoMPoST QCD KineticFree Streaming -- C -1 ∞ w~ Navier-Stokes Hydro -- 1-2/3 π w~ 0 0.2 0.4 0.6 0.8 1 1.2 0.01 0.1 1 10 E ff. s hea r- v i sc o s i t y : ( η / s ) e ff / ( η / s ) = / τ / τ R ( p T - p L ) / ( e + p ) Evolution time: w~=T( τ ) τ /(4 π η /s)Boltzmann RTA (T τ R =const)KoMPoST QCD KineticNavier-Stokes Hydro FIG. 2. (top) Evolution of energy-momentum tensor in therelaxation time approximation with constant η/s , comparedto Yang-Mills kinetic theory (KøMPøST) [44, 45]. (bottom)Evolution of the effective viscosity.
D. Evolution of energy momentum tensor forconformal system ( T ( τ ) τ R = const ) We now investigate the more commonly studied casewhere the relaxation time is chosen inversely proportionalto the eff. temperature, i.e. τ R T ( τ ) = 5( η/s ) = const ,as is appropriate for a conformal system [26, 35, 36, 39,56]. We conveniently introduce the dimensionless time-like variable x = ττ R = T ( τ ) τ η/s = 4 π w . (44)which one can also identify as the inverse Knudsen num-ber x ≡ Kn − . Since for the conformal system, the re-laxation time depends on the temperature of the systemas τ R = 5 η/sT ( τ ), it is convenient to perform a changeof variables to express τ ∂ τ = a ( x ) x∂ x , (45) where we defined the scale-factor a ( x ) = (cid:20)
23 + 14 (cid:18) b , + b , +2 C ( x ) C ( x ) (cid:19)(cid:21) , (46)and used the equation of motion for the energy densityalong with the fact that for an ultra-relativistic system dTT = de e . The evolution equation for the moments canbe then re-cast into the form a ( x ) x∂ x C ml = b ml, − C ml − + b ml, C ml + b ml, +2 C ml +2 − x (cid:16) C ml − C ml | eq (cid:17) . (47)indicating that for a given initial condition the non-equilibrium evolution of the system is uniquely deter-mined by the conformal scaling variable x .Numerical results for the evolution of the backgroundenergy-momentum tensor are presented in Fig. 2, wherewe also compare to the results for Yang-Mills kinetic the-ory obtained in [44, 45]. Clearly the overall behaviorof the curves is quite similar, showing a smooth tran-sition from an approximate free-streaming behavior atearly times towards the universal hydrodynamic behav-ior τ / e ( τ ) = (cid:16) τ / e (cid:17) ∞ (cid:18) − η/sT ( τ ) τ + · · · (cid:19) . (48)at late times. Similarly, to our previous discussion, thisbehavior can also be analyzed in terms of an effectiveshear viscosity over entropy ratio in far from equilibriumregime, which for the conformal system is given by( η/s ) eff η/s = 316 (cid:12)(cid:12)(cid:12)(cid:12) Re − Kn (cid:12)(cid:12)(cid:12)(cid:12) , (49)which is again constructed such that ( η/s ) eff η/s approachesunity in the limit τ /τ R (cid:29)
1. Even though this ratiodepicted in the lower panel of Fig. 2 magnifies the differ-ences between the results for Yang-Mills kinetic theory[44, 45] and the relaxation time approximation, the over-all differences are at most at the 10% level for interme-diate times ˜ w ∼ e ( τ ) ∝ /τ for x (cid:28) τ / e (cid:0) τ / e (cid:1) ∞ x (cid:28) = 1 C ∞ (cid:18) T ( τ ) τ πη/s (cid:19) / . (50)By inverting this relation, one then obtains the energydensity at late times as [46] (cid:16) τ / e (cid:17) ∞ = C ∞ (cid:18) πη/sT ( τ ) τ (cid:19) / ( eτ ) . (51)Specifically, for the conformal relaxation time approxi-mation, we find C ∞ ≈ . C ∞ ≈ IV. ENERGY MOMENTUM PERTURBATIONSAROUND BJORKEN FLOW
So far we have addressed the non-equilibrium evolutionof the average boost invariant and homogenous back-ground. We will now consider the propagation of lin- earized perturbations, sourced by (small) space-time de-pendent deviations of the initial energy momentum ten-sor from its (local) average. By linearizing the kineticequation around the boost invariant and homogenousbackground one finds an evolution equation for the per-turbation of the distribution function δf , i.e., (cid:104) p τ ∂ τ + p i ∂ i − p ς τ ∂ ς (cid:105) δf ( x, p ) = − p τ τ R δf ( x, p ) − p µ δu µ ( x ) τ R (cid:20)(cid:18) f BG ( τ, p T , | p ς | ) − f eq (cid:16) p τ T ( τ ) (cid:17)(cid:19) − p τ T ( τ ) f (cid:48) eq (cid:16) p τ T ( τ ) (cid:17)(cid:21) − p τ τ R δT ( x ) T ( τ ) (cid:20) p τ T ( τ ) f (cid:48) eq (cid:16) p τ T ( τ ) (cid:17) − T ( τ ) τ R ∂τ R ∂T (cid:18) f BG ( τ, p T , | p ς | ) − f eq (cid:16) p τ T ( τ ) (cid:17)(cid:19)(cid:21) . (52)where f (cid:48) eq ( x ) = df eq ( x ) /dx denotes the derivative of theequilibrium distribution. In the above expression, thechange in the rest-frame velocity δu µ ( x ) and local equi-librium temperature δT ( x ) are to be determined self-consistently from the (linearized) Landau matching con-dition. Starting from the linearized perturbations of theenergy momentum tensor δT µν ( x ) = (cid:104) p µ p µ (cid:105) δf (53)the change in the rest-frame velocity δu µ and energy den-sity in the local rest-frame δe are determined from thelinearized eigenvalue equation u µ δT µν + δu µ T µν = δeu ν + eδu ν , (54)with u µ δu µ = 0 . (55)By using the leading order solution u µ T µν = eu ν of theeigenvalue problem, one reads off δe = u µ δT µν u ν = δT ττ , (56a) δu τ = 0 , δu i = δT τi e + p T , δu ς = δT τς e + p L . (56b)Our strategy to determine the evolution of energy andmomentum perturbations then consist in determining allthe coefficients on the r.h.s. of Eq. (52) based on certainmoments of the perturbation of the phase-space distribu-tion function. In the following section we shall explainthis procedure in detail. A. Evolution equations for energy-momentumperturbations in the transverse plane
We will from now on restrict our attention to energy-momentum pertubations in the transverse plane, i.e. we will only consider variations in the transverse coordinates x . Since we are considering the evolution of linearizedperturbations on top of a (transversely) homogenousbackground, it is natural to express them in a Fourierbasis such that for each k -mode δf ( τ, x , p , | p ς | ) = (cid:90) d k (2 π ) δf k ( τ, p , | p ς | ) e i k · x , (57)where we denote δf k ( τ, p , | p ς | ) ≡ δf ( τ, k , p , | p ς | ). In ourapproach the Z subgroup symmetry of reflections alongthe beam axis is not broken at the level of the perturba-tions so the longitudinal rest frame velocity δu ς vanishesidentically. The transverse flow velocity can be decom-posed in the components parallel ( δu (cid:107) k ( τ )) and perpen-dicular ( δu ⊥ k ( τ )) to the wave vector k . Hence the pertur-bations to the thermodynamic fields take the followingform δT ( τ, x ) = (cid:90) d k (2 π ) δT k ( τ ) e i k · x , (58a) δu i ( τ, x ) = (cid:90) d k (2 π ) k j | k | (cid:104) δu (cid:107) k ( τ ) δ ji + δu ⊥ k ( τ ) (cid:15) ji (cid:105) e i k · x , (58b) δu τ ( x ) = 0 , δu ς ( x ) = 0 , (58c)where for physical perturbations, the reality conditions ofthe perturbations imply δT k ( τ ) = δT ∗− k ( τ ) and δu i k ( τ ) = δu i, ∗− k ( τ ). Denoting k · p | k | p τ = δ ij k i p j | k | p τ = cos( φ pk ) sin( θ p ) , (59) k × p | k | p τ = (cid:15) ij k i p j | k | p τ = sin( φ pk ) sin( θ p ) , (60)with φ pk = φ p − φ k being the angle between p and k in the transverse plane and sin θ p = | p | /p τ and insertingthe above expressions into the linearized kinetic equation(52) one then finds (cid:20) τ ∂ τ + i ( | k | τ ) k · p | k | p τ (cid:21) δf k ( τ, p , | p ς | ) = − (cid:18) ττ R (cid:19) δf k ( τ, p , | p ς | ) − (cid:18) ττ R (cid:19) (cid:20) δT k ( τ ) T ( τ ) + δu (cid:107) k ( τ ) k · p | k | p τ + δu ⊥ k ( τ ) k × p | k | p τ (cid:21) p τ T ( τ ) f (cid:48) eq (cid:16) p τ T ( τ ) (cid:17) + (cid:18) ττ R (cid:19) (cid:20) δu (cid:107) k ( τ ) k · p | k | p τ + δu ⊥ k ( τ ) k × p | k | p τ + δT k ( τ ) T ( τ ) T ( τ ) τ R ∂τ R ∂T (cid:21) (cid:18) f BG ( τ, | p | , | p ς | ) − f eq (cid:18) p τ T ( τ ) (cid:19)(cid:19) , (61)where the terms on the left hand side describe free-streaming, the terms in the first line correspond to therelaxation of the perturbation, the terms in the secondline describe the change of the equilibrium distributiondue to the perturbations and the terms in the last line de-scribe the change in the relaxation of the non-equilibriumbackground due to perturbations. B. Evolution equation of the spherical harmonicmoments
We follow the same strategy as for the evolution ofthe background and transform the evolution equation forthe distribution function into a coupled set of evolutionequations for the spherical harmonic moments δC ml, k ( τ ) = ν eff (cid:90) dp ς (2 π ) (cid:90) d p (2 π ) τ / p τ Y ml ( φ pk , θ p ) δf k ( τ, p , | p ς | ) , (62)where the azimuthal and polar angles φ pk , θ p are mea-sured with respect to the transverse wave vector ( k )and the longitudinal rapidity ( ς ) axis. Similarly as forthe background, the various components of the energy-momentum tensor are explicitly given in terms of thelowest order ( (cid:96) = 0 , ,
2) moments, as τ / δT ττ k = √ πδC , k , (63a) δ ij τ / δT ij k = (cid:114) π δC , k − (cid:114) π δC , k , (63e) k i k j k τ / δT ij k = (cid:114) π δC , k − (cid:114) π δC , k + (cid:114) π (cid:16) δC +22 , k + δC − , k (cid:17) , (63f) (cid:15) lj k i k l k τ / δT ij k = − i (cid:114) π (cid:16) δC +22 , k − δC − , k (cid:17) , (63g) δ ij i k i | k | τ / ( − τ ) δT ςj k = − i (cid:114) π (cid:16) δC +12 , k − δC − , k (cid:17) , (63h) (cid:15) ij i k i | k | τ / ( − τ ) δT ςj k = − (cid:114) π (cid:16) δC +12 , k + δC − , k (cid:17) , (63i) τ / τ δT ςς k = (cid:114) π δC , k + (cid:114) π δC , k , (63j)where we have conveniently decomposed all transversevectors δT τi , δT ςi into components parallel, perpendicu-lar and independent w.r.t. to the wave-vector k / | k | (andsimilarly for the tensor δT ij ). Due to the residual Z symmetry of long. reflection in rapidity, the components T τς and T ςi vanish identically, whereas all other compo-nents can in principle be non-zero in our setup.Evaluating the various couplings between sphericalharmonics using the identities listed in App. A, the evolu-tion equations for the spherical harmonic moments δC ml then take the form δ ij i k i | k | τ / δT τj k = − i (cid:114) π (cid:16) δC +11 , k − δC − , k (cid:17) , (63b) (cid:15) ij i k i | k | τ / δT τj k = − (cid:114) π (cid:16) δC +11 , k + δC − , k (cid:17) , (63c) τ / ( − τ ) δT τς k = (cid:114) π δC , k , (63d) Since physical perturbations of the phase-space distributions δf ( τ, x , p , | p ς | ) are real-valued, the linearized perturbations inFourier should also satisfy the condition δf − k ( τ, p , | p ς | ) = δf ∗ k ( τ, p , | p ς | ), such that in terms of the spherical harmonic mo-ments one finds δC m(cid:96), − k ( τ ) = ( − m (cid:16) δC − m(cid:96), k ( τ ) (cid:17) ∗ for physicalperturbations. τ ∂ τ δC ml, k = b ml, − δC ml − , k + b ml, δC ml, k + b ml, +2 δC ml +2 , k − i | k | τ (cid:16) u ml, − δC m +1 l − , k + u ml, + δC m +1 l +1 , k + d ml, − δC m − l − , k + d ml, + δC m − l +1 , k (cid:17) − (cid:18) ττ R (cid:19) (cid:20) δC ml, k + δe k e ( C (cid:48) eq ) ml (cid:21) − δe k e (cid:18) ττ R (cid:19) T ( τ ) τ R ∂τ R ∂T ( C eq − C ) ml − (cid:18) ττ R (cid:19) δu (cid:107) k (cid:16) u ml, − ( C eq − C + C (cid:48) eq ) m +1 l − + u ml, + ( C eq − C + C (cid:48) eq ) m +1 l +1 + d ml, − ( C eq − C + C (cid:48) eq ) m − l − + d ml, + ( C eq − C + C (cid:48) eq ) m − l +1 (cid:17) − (cid:18) ττ R (cid:19) δu ⊥ k i (cid:16) u ml, − ( C eq − C + C (cid:48) eq ) m +1 l − + u ml, + ( C eq − C + C (cid:48) eq ) m +1 l +1 − d ml, − ( C eq − C + C (cid:48) eq ) m − l − − d ml, + ( C eq − C + C (cid:48) eq ) m − l +1 (cid:17) (64)where the coefficients b ml are the same as in Eq. (28) and u ml, − = + (cid:114) ( l − m )( l − m − l + 1 , u ml, + = − (cid:115) ( l + m + 1)( l + m + 2)3 + 4 l ( l + 2) , (65) d ml, − = − (cid:114) ( l + m )( l + m − l + 1 , d ml, + = + (cid:115) ( l − m + 1)( l − m + 2)3 + 4 l ( l + 2) . (66)Physically the terms in the first line of Eq. (64) corre-spond to a free-streaming evolution, while the terms inthe last few lines capture the relaxation towards equilib-rium, including the changes of the equilibrium distribu-tion and background equilibration. By C (cid:48) eq in Eq. (64)we denote the moments( C (cid:48) eq ) ml ( τ ) = (cid:90) dp ς (2 π ) (cid:90) d p (2 π ) τ / p τ (67) Y ml ( φ pk , θ p ) (cid:18) p τ T ( τ ) (cid:19) f (cid:48) eq (cid:18) p τ T ( τ ) (cid:19) , which can be determined via integration by parts as( C (cid:48) eq ) ml ( τ ) = − C eq ) ml , (68)with the equilibrium moments ( C eq ) ml determined byEq. (26). By inserting the relations in Eqns. (63) intothe linearized Landau matching conditions in Eq. (56),one also finds that the linearized perturbations of the en-ergy density δe k and long. and transverse flow velocities δu (cid:107) k and δu ⊥ k are given by τ / δe k ( τ ) = √ πδC , k ( τ ) , (69) τ / ( e + p T ) δu (cid:107) k ( τ ) = − (cid:114) π (cid:16) δC , k ( τ ) − δC − , k ( τ ) (cid:17) ,τ / ( e + p T ) δu ⊥ k ( τ ) = i (cid:114) π (cid:16) δC , k ( τ ) + δC − , k ( τ ) (cid:17) , such that the evolution equations for the moments inEq. (64) form a closed set of equations. We further note that by virtue of the decomposition into spherical har-monic moments, Y ml ( φ pk , θ p ), the information on the di-rection of the transverse wave-vector k has disappearedfrom the evolution equation, which as a consequence ofthe azimuthal rotation symmetry of the background onlydepends on the magnitude of the wave-vector | k | .Before we address the physically relevant initial condi-tions for the moments δC ml, k , we also note that it is oftenuseful to consider the evolution equation at a fixed valueof propagation phase κ = | k | ( τ − τ ) (70)rather than a fixed value for the wave-number | k | . Bychanging the variables from | k | to κ = | k | ( τ − τ ) forthe moments, the time derivate needs to be evaluatedaccording to τ ∂ τ | k = τ ∂ τ | k ( τ − τ ) + ττ − τ | k | ( τ − τ ) ∂ | k | ( τ − τ ) (cid:12)(cid:12) τ , (71)such that the evolution equation for the moments re-ceives one additional term associated with this changeof variables. Similarly, for a conformal system, it is alsoconvenient to express the evolution in terms of the scaledevolution time x = τ /τ R introduced in Eq. (44). Startingfrom Eq. (64) and denoting s ( τ ) = τ − τ τ , with a ( x ) x∂ x s ( x ) = (1 − s ( x )) , (72)to perform these changes, the equation of motion for themoments then takes the form0 (cid:104) s ( x ) a ( x ) x∂ x + κ∂ κ (cid:105) δC ml,κ = s ( x ) (cid:16) b ml, − δC ml − ,κ + b ml, δC ml,κ + b ml, +2 δC ml +2 ,κ (cid:17) − iκ (cid:16) u ml, − δC m +1 l − ,κ + u ml, + δC m +1 l +1 ,κ + d ml, − δC m − l − ,κ + d ml, + δC m − l +1 ,κ (cid:17) − xs ( x ) (cid:20) δC ml,κ + δe κ e ( C (cid:48) eq ) ml (cid:21) − xs ( x ) δe κ e T ( τ ) τ R ∂τ R ∂T ( C eq − C ) ml − xs ( x ) δu (cid:107) κ (cid:16) u ml, − ( C eq − C + C (cid:48) eq ) m +1 l − + u ml, + ( C eq − C + C (cid:48) eq ) m +1 l +1 + d ml, − ( C eq − C + C (cid:48) eq ) m − l − + d ml, + ( C eq − C + C (cid:48) eq ) m − l +1 (cid:17) − xs ( x ) δu ⊥ κ i (cid:16) u ml, − ( C eq − C + C (cid:48) eq ) m +1 l − + u ml, + ( C eq − C + C (cid:48) eq ) m +1 l +1 − d ml, − ( C eq − C + C (cid:48) eq ) m − l − − d ml, + ( C eq − C + C (cid:48) eq ) m − l +1 (cid:17) (73)which we will employ below to obtain the numerical so-lution for the evolution of perturbations.Based on Eq. (73), one also explicitly observes thatin the limit τ /τ R → s ( x ) = 1 is a fixed pointof Eq. (72), the solution to the evolution equation forthe moments δC ml, k only depends on the propagationphase κ = | k | ( τ − τ ) and the scaled evolution time x = τ /τ R = T ( τ ) τ η/s . While this conformal scaling behav-ior was empirically observed in [44, 45] from numericalsolutions of the Boltzmann equation in QCD kinetic the-ory at different values of the coupling strength λ = g N c ,it is interesting to point out that in the present contextthe conformal scaling behavior directly manifests itself atthe level of the equations of motion. C. Initial conditions for energy and momentumperturbations
So far we have discussed, the evolution equations forlinearized perturbations on top of a Bjorken background.Now in order to apply this framework to describe theearly time dynamics of high-energy heavy-ion collisions,the equations of motion need to be supplemented bysuitable initial conditions, which describe the associatedchange of the phase-space density at initial time. Whilein principle one could imagine a large variety of differentinitial conditions, we will follow [44, 45] and only considerthe response of the system to changes of the conservedquantities of the system, associated with initial energyand momentum perturbations as detailed below.
1. Energy perturbations
We follow the arguments of [44] and associate initialenergy perturbations with an infinitesimal change of theenergy scale of the background distribution, such thatthe associated phase-space distribution for energy per- turbations is given by δf k ( τ , p , | p ς | ) = − (cid:18) | p | ∂ | p | f (0) BG ( | p | , p ς ) (cid:19) e − i k p | p | τ , (74)where we introduced the phase-factor e − i k p | p | τ to ac-count for the free-streaming evolution at early times τ < τ (cid:28) τ R . Based on the explicit form of the ini-tial background distribution in Eq. (29), the integrals forthe moments δC ml in Eq. (62) can be evaluated using12 π (cid:90) π dφ pk e − i | k | τ cos( φ pk ) e imφ pk = ( − i ) m J m ( | k | τ ) , (75)along with Eq. (31) yielding δC ml, k ( τ ) = τ / ( eτ ) ( − i ) m J m ( | k | τ ) y ml P ml (0) , (76)where ( eτ ) denotes the asymptotic energy density of thebackground (c.f. Sec. III). Specifically, for the energy andmomentum density one has δe k ( τ ) e = J ( | k | τ ) , (77)( e + p T ) e δu (cid:107) k ( τ ) = − iJ ( | k | τ ) , (78)( e + p T ) e δu ⊥ k ( τ ) = 0 , (79)reproducing the result of [44] for the free-streaming re-sponse function.
2. Momentum perturbations
Similarly, we associate initial momentum perturba-tions with an infinitesimal change of the transverse ve-locity of the background, such that following the argu-ments of [44] the associated phase-space distribution formomentum perturbations is given by δf i k ( τ , p , p ς ) = 2 (cid:18) p i ∂ | p | f (0) BG ( | p | , p ς ) (cid:19) e − i k p | p | τ , (80)1where the index i contains the information about thedirection of the initial momentum perturbation. Decom-posing δf i k into the directions parallel and perpendicularto the wave-vector, we can distinguish between longitu-dinal and transverse momentum perturbations δf (cid:107) k ( τ , p , p ς ) = (81)2 cos( φ pk ) (cid:18) | p | ∂ | p | f (0) BG ( | p | , p ς ) (cid:19) e − i k p | p | τ ,δf ⊥ k ( τ , p , p ς ) = (82)2 sin( φ pk ) (cid:18) | p | ∂ | p | f (0) BG ( | p | , p ς ) (cid:19) e − i k p | p | τ . Evaluating the moments, one then finds δC m (cid:107) l, k ( τ ) = (83) − iτ / ( eτ ) ( − i ) m (cid:16) J m +1 ( | k | τ ) − J m − ( | k | τ ) (cid:17) y ml P ml (0) , for longitudinal momentum perturbations and similarlyfor transverse momentum perturbations δC m ⊥ l, k ( τ ) = (84) − τ / ( eτ ) ( − i ) m (cid:16) J m +1 ( | k | τ ) + J m − ( | k | τ ) (cid:17) y ml P ml (0) . Specifically, for the energy and momentum density onehas δe k ( τ ) e = − iJ ( | k | τ ) , (85)( e + p T ) e δu (cid:107) k ( τ ) = 2 J ( | k | τ ) | k | τ − J ( | k | τ ) , (86)( e + p T ) e δu ⊥ k ( τ ) = 0 , (87)for longitudinal momentum perturbations, whereas fortransverse momentum perturbations δe k ( τ ) e = 0 , (88)( e + p T ) e δu (cid:107) k ( τ ) = 0 , (89)( e + p T ) e δu ⊥ k ( τ ) = 2 J ( | k | τ ) | k | τ , (90)in agreement with [44]. V. NON-EQUILIBRIUM GREEN’S FUNCTIONSOF THE ENERGY-MOMENTUM TENSOR
We now proceed to the calculation of the response ofthe energy-moment tensor to initial energy and momen-tum perturbations, based on numerical solutions of theevolution equations for the moments δC ml, k starting fromthe initial condition described in the previous section.Since the set of moments δC ml, k contain an overwhelmingamount of information, we will therefore restrict our at-tention to the evolution of the low order moments δC ml, k with (cid:96) ≤
2, which can be related to the various compo-nents of the energy momentum tensor δT µν k according toEq. (63). Instead of investigating the dynamics of indi-vidual moments δC ml, k directly, we find it more insightfulto consider the linear response functions G µναβ ( k , τ, τ ) in-troduced in [44, 45], such that δT µν k ( τ ) e ( τ ) = 12 G µναβ ( k , τ, τ ) δT αβ k ( τ ) e ( τ ) . (91)Noteably, the Green’s functions G µναβ ( k , τ, τ ) provide thebuiliding block of the pre-equilibrium computer codeKøMPøST [63] and can be obtained in terms of linearcombinations of the moments as described below. Sincewe are primarily interested in the limit τ /τ R →
0, wherethe kinetic theory framework describes the equilibrationprocess of the system all the way from very early timesup to the onset of hydrodynamic behavior, we will dropthe explicit dependence on τ in the following to lightenthe notation. By expressing the response to an initialenergy perturbation in the following basis of scalars (s),vectors (v) and tensors (t) G ττττ ( k , τ ) = G ss ( κ, x ) , (92a) G τiττ ( k , τ ) = − i k i | k | G vs ( κ, x ) , (92b) G ijττ ( k , τ ) = δ ij G t,δs ( κ, x ) + k i k j k G t,ks ( κ, x ) , (92c)and adapting the normalization δe ( τ ) /e ( τ ) = 1, therelevant response functions can then be determined from(c.f. Appendix B of [44]) G ss ( κ, x ) = δe κ ( x ) e ( x ) , (93a) G vs ( κ, x ) = δ ij i k i | k | δT jκ ( x ) e ( x ) , (93b) G t,δs ( κ, x ) = (cid:20) δ ij − k i k j k (cid:21) δT ijκ ( x ) e ( x ) , (93c) G t,ks ( κ, x ) = (cid:20) k i k j k − δ ij (cid:21) δT ijκ ( x ) e ( x ) . Similarly, the response to an initial momentum pertur-bation can be characterized by a set of six independentresponse functions, G τττl ( k , τ ) = − i k l | k | G sv ( κ, x ) , (94a) G τiτl ( k , τ ) = δ ij G δv ( κ, x ) + k i k j k G kv ( κ, x ) , (94b)2 -1.5-1-0.5 0 0.5 1 1.5 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) E ne r g y r e s pon s e : G ss ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12-0.6-0.4-0.2 0 0.2 0.4 0.6 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) M o m en t u m r e s pon s e : G sv ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12 -1.5-1-0.5 0 0.5 1 1.5 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) P r e ss u r e r e s pon s e : G s t, δ ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12-0.6-0.4-0.2 0 0.2 0.4 0.6 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) S hea r- s t r e ss r e s pon s e : G s t, k ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12 FIG. 3. Evolution of spectrum of energy-momentum perturbations in response to an initial energy perturbation. Differentcurves in each panel correspond to different evolution times T ( τ ) τ / (4 πη/s ); different panels show the response of the differentcomponents of the energy-momentum tensor as a function of the wave-number | k | ( τ − τ ). G ijτl ( k , τ ) = − i k l | k | δ ij G t,δv ( κ, x ) − i k i k j k k | k | G t,kv ( κ, x ) − i δ il k j + δ jl k i | k | G t,mv ( κ, x ) (94c)which upon adapting the normalization condition δT τiτj,κ /e ( τ ) = δ ij , where as in Eq. (91) upper indices( τ i ) refer to the respective components of the energy-momentum tensor and lower indices ( τ j ) indicate the di-rection of the initial momentum perturbation, are givenby (c.f. Appendix B of [44]) G sv ( κ, x ) = δ ij i k i | k | δT τττj,κ ( x ) δe ( x ) , (95a) G v,δv ( κ, x ) = (cid:20) δ ij − k i k j k (cid:21) δT τiτj,κ ( x ) δe ( x ) , (95b) G v,kv ( κ, x ) = (cid:20) k i k j k − δ ij (cid:21) δT τiτj,κ ( x ) δe ( x ) , (95c) G t,δv ( κ, x ) = (cid:20) δ ij − k i k j k (cid:21) i k l | k | δT ijτl,κ ( x ) δe ( x ) (95d) G t,mv ( κ, x ) = 2 i k i | k | (cid:20) δ jl − k j k l k (cid:21) δT ijτl,κ ( x ) δe ( x ) (95e) G t,kv ( κ, x ) = (95f) (cid:18)(cid:20) k i k j k − δ ij (cid:21) i k l | k | − i k i | k | (cid:20) δ jl − k j k l k (cid:21)(cid:19) δT ijτl,κ ( x ) δe ( x ) . While Eqns. (92) and (94) provide a basis for ex-pressing the response of the energy-momentum tensor G µναβ ( k , τ, τ ) in wave-number ( k ) space, we will also beinterested in the response of the energy-momentum ten-sor in coordinate space G µναβ ( x − x , τ, τ ), where an analo-gous decomposition can be performed w.r.t. to the vector x − x . We will refrain from presenting all the details andinstead refer the interested reader to Appendix B of [44].3 -4-3-2-1 0 1 2 3 4 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) E ne r g y r e s pon s e : G vs ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12-4-3-2-1 0 1 2 3 4 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) M o m en t u m r e s pon s e : G vv , δ ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12-1-0.5 0 0.5 1 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) S hea r- s t r e ss r e s pon s e : G v t, m ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12 -4-3-2-1 0 1 2 3 4 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) P r e ss u r e r e s pon s e : G v t, δ ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12-4-3-2-1 0 1 2 3 4 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) M o m en t u m r e s pon s e : G vv , k ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12-1-0.5 0 0.5 1 0 5 10 15 20 25 30 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) S hea r- s t r e ss r e s pon s e : G v t, k ( w ~ , k ∆ τ ) Wave number: k ∆τ free streaming 12 FIG. 4. Evolution of spectrum of energy-momentum perturbations in response to an momentum energy perturbation. Differentcurves in each panel correspond to different evolution times T ( τ ) τ / (4 πη/s ); different panels show the response of the differentcomponents of the energy-momentum tensor as a function of the wave-number | k | ( τ − τ ). A. Numerical results
We will now present numerical results for the variousresponse functions, focusing on the case of a conformalrelaxation time. We follow the same strategy as for thebackground and solve a truncated set of evolution equa- tions, checking that including higher order moments doesnot significantly alter the results. Noteably, we find thata rather large set of moments is required to achieve appar-ent convergence and we will present results for l max = 64in the following.Our results for the evolution of the response func-4 -0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) E ne r g y r e s pon s e : ∆ τ G ss ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) M o m en t u m r e s pon s e : ∆ τ G sv ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012 -0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) P r e ss u r e r e s pon s e : ∆ τ G s t, δ ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) S hea r- s t r e ss r e s pon s e : ∆ τ G s t, x ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012 FIG. 5. Evolution of the energy-momentum response to an initial energy perturbation in coordinate space, based on RTA (solid)and Yang-Mills kinetic theory (KøMPøST)(dotted) [44] . Different curves in each panel correspond to different evolution times T ( τ ) τ / (4 πη/s ); different panels show the response of the different components of the energy-momentum tensor as a functionof the propagation distance | x − x | / ( τ − τ ). tions are compactly summarized in Figs. 3 and 4 , wherewe present the spectrum of perturbations G µνττ ( k , τ − τ , τ /τ R ) and G µντi ( k , τ − τ , τ /τ R ) as a function of | k | ( τ − τ ) for different values of the evolution time˜ w = T ( τ ) τ / (4 πη/s ). Different panels in Figs. 3 and 4show the results for the different ( µν ) components, de-composed in the tensor basis described above. In orderto facilitate the interpretation, we have also labeled thevarious response functions according to the componentsof the energy-momentum tensor that they affect.Starting from the initial free-streaming behavior atearly times T ( τ ) τ / (4 πη/s ) (cid:28) | k | ( τ − τ ) modes sets in, such that bythe time the system enters the hydrodynamic regime( ˜ w = T ( τ ) τ / (4 πη/s ) ∼
1) only the long wave-lengthmodes survive. Since at early times the system is highlyanisotropic, the longitudinal pressure is effectively zeroand transverse perturbations initially propagate with aphase-velocity of nearly the speed of light. Subsequently, as the system becomes more and more isotropic thephase-velocity decreases and eventually approaching thespeed of sound, which in Figs. 3 and 4 results in a shift ofthe peaks towards larger values of the propagation phase | k | ( τ − τ ). Strikingly, the qualitative behavior observedfrom Figs. 3 and 4 is very similar to the results obtainedin Yang-Mills kinetic theory in [44], albeit we find thatin the relaxation time approximation the viscous damp-ing of short wave length modes becomes efficient on asomewhat shorter time scale.Even though some of the features of the evolution canbe understood quite naturally in wave-number ( k ) space,in practice one is mostly interested in the Green’s func-tions G µνττ ( x − x , τ − τ , τ /τ R ) in position space, whichdirectly describe the physical response of the energy mo-mentum tensor δT µν ( τ, x ) to a localized initial energyperturbation δT ττ ( τ , x ) according to δT µν ( τ, x ) e ( τ ) = 12 (cid:90) x G µναβ ( x − x , τ − τ ) δT αβ ( τ , x ) e ( τ ) . (96)5 -2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) E ne r g y r e s pon s e : ∆ τ G vs ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012-2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) M o m en t u m r e s pon s e : ∆ τ G vv , δ ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012-2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) S hea r- s t r e ss r e s pon s e : ∆ τ G s t, m ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012 -2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) P r e ss u r e r e s pon s e : ∆ τ G v t, δ ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012-2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) M o m en t u m r e s pon s e : ∆ τ G vv , x ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012-2-1.5-1-0.5 0 0.5 1 1.5 2 2.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 E v o l u t i on t i m e : w ~= T ( τ ) τ / ( π η / s ) S hea r- s t r e ss r e s pon s e : ∆ τ G s t, x ( ∆ x , ∆ τ , w ~ ) Distance: ∆ x/ ∆τ free streamingw~=0.25w~=0.50w~=1.00w~=1.50w~=2.00 012 FIG. 6. Evolution of the energy-momentum response to an initial momentum perturbation in coordinate space, based onRTA (solid) and Yang-Mills kinetic theory (KøMPøST)(dotted) [44]. Different curves in each panel correspond to differentevolution times T ( τ ) τ / (4 πη/s ); different panels show the response of the different components of the energy-momentum tensoras a function of the propagation distance | x − x | / ( τ − τ ). In practice, the coordinate space response can be ob-tained in a straightforward way via a set of Bessel-Fourier transforms of the response functions, and weuse the KøMPøSTsoftware [63] to perform this task.We note that when implementing our results into KøMPøST, one needs to take into account that the re-sults in [44] are presented in terms of the time vari-able x Id S = T Id ( τ ) τ / (4 πη/s ), where T Id ( τ ) is definedvia the asymptotic temperature according to T Id ( τ ) = τ − / lim τ →∞ ( T τ / ). However, from the point of view6of the RTA Boltzmann equation it is more natural tostudy the evolution as a function of the scaling variable˜ w = T ( τ ) τ / (4 πη/s ), where T ( τ ) denotes the equilib-rium temperature obtained from e ( τ ) via Landau match-ing. Of course, the two quantities are related by x Id S =˜ w (cid:16) ( τ / e ) ∞ τ / e ( τ ) (cid:17) / and we have taken this difference intoaccount in all explicit comparisons presented in this pa-per. In order to provide an apples to apples comparisonbetween results obtained within the relaxation time ap-proximation and Yang-Mills kinetic theory, the differentresponse functions are smeared out with the same smear-ing kernel K σ = exp( − k ( τ − τ ) / σ ) as in [44], whichfor Yang-Mills kinetic theory results is necessary in orderto stabilize the numerical Bessel-Fourier transform.Our results for the coordinate space response functionsare shown in Figs. 5 and 6, where we display the variousresponse functions as a function of the propagation dis-tance | x − x | / ( τ − τ ). Solid curves in Figs. 5 and 6 showour results obtained from the Boltzmann equation in re-laxation time approximation, which are compared to theresults obtained in Yang-Mills kinetic theory from [44]shown as dashed curves. Some of the most important fea-tures that can be immediately observed from Figs. 5 and6, include the viscous broadening of the peaks due to thedamping of high wave-number modes, as well as the shiftof the peaks towards smaller values of | x − x | / ( τ − τ )associated with the aforementioned change of the effec-tive (transverse) speed of sound, due to the increase ofthe longitudinal pressure. One also observes a decreaseof the shear-stress response compared to the pressure re-sponse (for both energy and momentum perturbations),such that by ˜ w = T ( τ ) τ / (4 πη/s ) ∼ ∼
10% on the relevant time scales˜ w = T ( τ ) τ / (4 πη/s ) ∼ VI. CONCLUSIONS & OUTLOOK
We derived a new method to calculate non-equilibriumGreen’s function of the energy momentum tensor basedon moment equations of kinetic equations for linearizedperturbations. Due to the particularly simple structureof the relaxation time approximation considered in thiswork, we obtained a closed set of moment equations for the evolution of the dimension four moments C ml and δC ml , which are relevant to study the evolution of theenergy-momentum tensor, that is determined by the low-est order ( (cid:96) ≤
2) moments. Even though, for morecomplex interactions, a truncation at the level of dimen-sion four operators is no longer sufficient to obtain aclosed set of evolution equations, we naturally believethat by including higher order operators our method canbe extended to systems with more complex interactions,thereby generalizing the usual moment method at levelof perturbations.Based on the evolution equations for the moments C ml and δC ml in Eqns. (27) and (64), we studied the evo-lution of average energy-momentum tensor in Bjorkenflow, as well as the out-of-equilibrium linear response ofthe system to initial energy and momentum perturba-tions in the transverse plane. By truncating the infi-nite hierarchy of moment equations at large finite order,we obtained numerical solutions for the evolution of thenon-equilibrium background and the Green’s functions.When comparing our results to previous calculations ofKøMPøSTin Yang-Mills theory [44], we found a strikingsimilarity between the different theories. Even thoughthe macroscopic differences between the two microscopiccalculations are only at the ten percent level, it wouldbe interesting to explore to what extent these can affectconcrete observables, such as e.g. the flow harmonics v n ,the charged particle multiplicity dN ch /dη or the trans-verse energy dE ⊥ /dη . Since there are first hints thatin particular the charged particle multiplicity dN ch /dη ,may be a rather sensitive measure of the entropy produc-tion during the pre-equilibrium phase [46], this remainsan interesting question which we expect to be addressedin more detail in future studies.Besides our numerical studies, we also found that var-ious conformal scaling features, which have been empir-ically observed in [44], can be directly seen at level ofthe equations of motion for the moments δC ml , and weexpect that further analytic insights into the structureof the non-equilibrium Green’s functions. Specifically, itwould be interesting to further explore for example theearly and late time asymptotics of the Green’s functionsbased on this formulation. Furthermore, one can also in-vestigate to what extent the highly non-trivial evolutionof the Green’s functions can be captured by ”renormal-ized transport coefficients” as suggested in various works[36, 53, 64] studying the highly symmetric Bjorken flow.Beyond such analytic insights, it would also interestingto further systematically extent the pre-equilibrium de-scription in terms of non-equilibrium Green’s functions.Since the numerical solution of the moment equations iscomparatively straightforward, the Boltzmann equationin relaxation time approximation provides an idealtesting ground for such ideas. Specifically, with themethodology developed in this paper it should be com-paratively straightforward to extent the description toinclude also longitudinal fluctuations, investigate Green’sfunctions for higher order moments or study the effect7of additional conserved charges, all of which are ratherchallenging tasks within a QCD kinetic description. Byexplicitly comparing the linearized description in termsof non-equilibrium Green’s functions to full numericalsolutions of the RTA Boltzmann equation, one couldalso obtain additional insights into the reliability andbreakdown of the linearized description and potentiallyimprove the range of applicablity to describe smallcollision systems. ACKNOWLEDGMENTS
We thank A. Behtash for his involvement in thisproject during its early stages. MM and SK were sup-ported in part by the US Department of Energy grantDE-FG02-03ER41260. MM is also supported by theBEST (Beam Energy Scan Theory) DOE Topical Collab-oration. SS, PP and SO are supported by the DeutscheForschungsgemeinschaft (DFG, German Research Foun-dation) through the CRC-TR 211 “Strong-interactionmatter under extreme conditions” under Project number315477589.
Appendix A: Identities for spherical harmonics &associated Legendre polynomials
Below we summarize some of the identities used to de-rive the evolution equations for the background moments C ml and the linearized perturbations δC ml . Specifically,for the evolution equations of the background momentswe make use of the identities(1 − x ) ddx P ml ( x ) = ∆ ml, − P ml − ( x ) + ∆ ml, + P ml +1 , (A1)and xP ml ( x ) = ξ ml, − P ml − ( x ) + ξ ml, + P ml +1 ( x ) , (A2)as well as x P ml ( x ) = ξ (2) ,ml, − P ml − ( x ) + ξ (2) ,ml, P ml ( x ) + ξ (2) ,ml, +2 P ml +2 ( x ) , (A3)to derive the relation (cid:20)(cid:18) − x (cid:19) − x (1 − x ) ddx (cid:21) P ml ( x ) = (A4) a ml, − P ml − ( x ) + a ml, P ml ( x ) + a ml, +2 P ml +2 ( x ) . such that τ ∂ τ C ml = b ml, − C ml − + b ml, C ml + b ml, +2 C ml +2 (A5)+ (cid:90) d p (2 π ) (cid:90) dp ς (2 π ) (cid:16) τ / p τ (cid:17) Y ml ( φ, θ ) τ ∂ τ f ( p ) . which is the identity used in the main text. Similarly,making use of the relations P − ml ( x ) = σ ml P ml ( x ) , σ ml = ( − m ( l − m )!( l + m )! , (A6)along with (cid:112) − x P ml = 12 l + 1 (cid:16) P m +1 l − ( x ) − P m +1 l +1 ( x ) (cid:17) , (A7)we can evaluate the additional termssin( θ ) e + iφ Y ml ( φ, θ ) = u ml, − Y m +1 l − ( φ, θ ) + u ml, + Y m +1 l +1 ( φ, θ ) , (A8)sin( θ ) e − iφ Y ml ( φ, θ ) = d ml, − Y m − l − ( φ, θ ) + d ml, + Y m − l +1 ( φ, θ ) , (A9)which arise in the evolution equations for the linearizedenergy-momentum perturbations. Below we list the co-efficients entering the above identities∆ ml, − = ( l + 1)( l + m )2 l + 1 , (A10)∆ ml, + = − l ( l − m + 1)2 l + 1 ,ξ ml, − = l + m l + 1 , (A11) ξ ml, + = l − m + 12 l + 1 ,ξ (2) ,ml, − = ξ ml, − ξ ml − , − , (A12) ξ (2) ,ml, = (cid:0) ξ ml, − ξ ml − , + + ξ ml, + ξ ml +1 , − (cid:1) ,ξ (2) ,ml, +2 = ξ ml, + ξ ml +1 , + ,a ml, − = − ξ (2) ,ml, − − ∆ ml, − ξ ml − , − , (A13) a ml, = 13 − ξ (2) ,ml, − ∆ ml, − ξ ml − , + − ∆ ml, + ξ ml +1 , − ,a ml, +2 = − ξ (2) ,ml, +2 − ∆ ml, + ξ ml +1 , + .b ml, − = a ml, − y ml y ml − , (A14) b ml, = a ml, ,b ml, +2 = a ml, +2 y ml y ml +2 .u ml, − = + y ml (2 l + 1) y l − ,m +1 , (A15) u ml, + = − y ml (2 l + 1) y l +1 ,m +1 , d ml, − = + y ml (2 l + 1) y l − ,m − σ ml σ − m +1 l − , (A16) d ml, = − y ml (2 l + 1) y l +1 ,m − σ ml σ − m +1 l +1 , which upon further simplifications yield the resultsquoted in the main text. [1] J. Adams et al. (STAR), Nucl. Phys. A757 , 102 (2005),arXiv:nucl-ex/0501009 [nucl-ex].[2] K. Adcox et al. (PHENIX), Nucl. Phys.
A757 , 184(2005), arXiv:nucl-ex/0410003 [nucl-ex].[3] B. B. Back et al. , Nucl. Phys.
A757 , 28 (2005),arXiv:nucl-ex/0410022 [nucl-ex].[4] I. Arsene et al. (BRAHMS), Nucl. Phys.
A757 , 1 (2005),arXiv:nucl-ex/0410020 [nucl-ex].[5] U. Heinz and R. Snellings, Ann. Rev. Nucl. Part. Sci. ,123 (2013), arXiv:1301.2826 [nucl-th].[6] M. Luzum and H. Petersen, J. Phys. G41 , 063102 (2014),arXiv:1312.5503 [nucl-th].[7] C. Gale, S. Jeon, and B. Schenke, Int. J. Mod. Phys.
A28 , 1340011 (2013), arXiv:1301.5893 [nucl-th].[8] R. Derradi de Souza, T. Koide, and T. Kodama, Prog.Part. Nucl. Phys. , 35 (2016), arXiv:1506.03863 [nucl-th].[9] D. A. Teaney, in Quark-gluon plasma 4 , edited byR. C. Hwa and X.-N. Wang (2010) pp. 207–266,arXiv:0905.2433 [nucl-th].[10] D. Kharzeev, E. Levin, and M. Nardi, Phys. Rev.
C71 ,054903 (2005), arXiv:hep-ph/0111315 [hep-ph].[11] J. S. Moreland, J. E. Bernhard, and S. A. Bass, Phys.Rev.
C92 , 011901 (2015), arXiv:1412.4708 [nucl-th].[12] B. Alver, M. Baker, C. Loizides, and P. Steinberg,(2008), arXiv:0805.4411 [nucl-ex].[13] A. Kurkela and E. Lu, Phys. Rev. Lett. , 182301(2014), arXiv:1405.6318 [hep-ph].[14] W. van der Schee, P. Romatschke, and S. Pratt, Phys.Rev. Lett. , 222302 (2013), arXiv:1307.2539 [nucl-th].[15] B. Schenke, P. Tribedy, and R. Venugopalan, Phys. Rev.Lett. , 252301 (2012), arXiv:1202.6646 [nucl-th].[16] J. Nagle and W. Zajc, Phys. Rev. C , 054908 (2019),arXiv:1808.01276 [nucl-th].[17] M. Martinez, M. D. Sievert, D. E. Wertepny, andJ. Noronha-Hostler, (2019), arXiv:1911.10272 [nucl-th].[18] M. Martinez, M. D. Sievert, D. E. Wertepny, andJ. Noronha-Hostler, (2019), arXiv:1911.12454 [nucl-th].[19] S. Schlichting and D. Teaney, Ann. Rev. Nucl. Part. Sci. , 447 (2019), arXiv:1908.02113 [nucl-th].[20] P. Romatschke, Phys. Rev. Lett. , 012301 (2018),arXiv:1704.08699 [hep-th].[21] W. Florkowski, M. P. Heller, and M. Spalinski, Rept.Prog. Phys. , 046001 (2018), arXiv:1707.02282 [hep-ph].[22] P. Romatschke and U. Romatschke, Relativistic FluidDynamics In and Out of Equilibrium , Cambridge Mono-graphs on Mathematical Physics (Cambridge UniversityPress, 2019) arXiv:1712.05815 [nucl-th].[23] A. Kurkela and Y. Zhu, Phys. Rev. Lett. , 182301(2015), arXiv:1506.06647 [hep-ph].[24] R. Critelli, R. Rougemont, and J. Noronha, JHEP ,029 (2017), arXiv:1709.03131 [hep-th].[25] G. S. Denicol, U. W. Heinz, M. Martinez, J. Noronha,and M. Strickland, Phys.Rev.Lett. , 202301 (2014), arXiv:1408.5646 [hep-ph].[26] W. Florkowski, R. Ryblewski, and M. Strickland,Nucl.Phys. A916 , 249 (2013), arXiv:1304.0665 [nucl-th].[27] W. Florkowski, R. Ryblewski, and M. Strickland, Phys.Rev.
C88 , 024903 (2013), arXiv:1305.7234 [nucl-th].[28] G. S. Denicol, U. W. Heinz, M. Martinez, J. Noronha,and M. Strickland, Phys.Rev.
D90 , 125026.[29] P. M. Chesler and L. G. Yaffe, Phys.Rev.
D82 , 026006(2010), arXiv:0906.4426 [hep-th].[30] M. P. Heller, R. A. Janik, and P. Witaszczyk,Phys.Rev.Lett. , 201602 (2012), arXiv:1103.3452[hep-th].[31] W. van der Schee, Phys.Rev.
D87 , 061901 (2013),arXiv:1211.2218 [hep-th].[32] P. M. Chesler, JHEP , 146 (2016), arXiv:1601.01583[hep-th].[33] M. Martinez and M. Strickland, Nucl. Phys. A848 , 183(2010), arXiv:1007.0889 [nucl-th].[34] W. Florkowski and R. Ryblewski, Phys.Rev.
C83 , 034907(2011), arXiv:1007.0130 [nucl-th].[35] J.-P. Blaizot and L. Yan, Annals Phys. , 167993(2020), arXiv:1904.08677 [nucl-th].[36] J.-P. Blaizot and L. Yan, Phys. Lett.
B780 , 283 (2018),arXiv:1712.03856 [nucl-th].[37] J.-P. Blaizot and L. Yan, JHEP , 161 (2017),arXiv:1703.10694 [nucl-th].[38] A. Kurkela, W. van der Schee, U. A. Wiedemann,and B. Wu, Phys. Rev. Lett. , 102301 (2020),arXiv:1907.08101 [hep-ph].[39] M. P. Heller, A. Kurkela, M. Spaliski, and V. Svensson,Phys. Rev. D97 , 091503 (2018), arXiv:1609.04803 [nucl-th].[40] M. P. Heller and M. Spalinski, Phys. Rev. Lett. ,072501 (2015), arXiv:1503.07514 [hep-th].[41] S. Jaiswal, C. Chattopadhyay, A. Jaiswal, S. Pal,and U. Heinz, Phys. Rev. C , 034901 (2019),arXiv:1907.07965 [nucl-th].[42] J. Casalderrey-Solana, N. I. Gushterov, and B. Meiring,JHEP , 042 (2018), arXiv:1712.02772 [hep-th].[43] L. Keegan, A. Kurkela, A. Mazeliauskas, and D. Teaney,JHEP , 171 (2016), arXiv:1605.04287 [hep-ph].[44] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlicht-ing, and D. Teaney, Phys. Rev. C99 , 034910 (2019),arXiv:1805.00961 [hep-ph].[45] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlicht-ing, and D. Teaney, Phys. Rev. Lett. , 122302 (2019),arXiv:1805.01604 [hep-ph].[46] G. Giacalone, A. Mazeliauskas, and S. Schlichting, Phys.Rev. Lett. , 262301 (2019), arXiv:1908.02866 [hep-ph].[47] C. Gale, J.-F. Paquet, B. Schenke, and C. Shen, in (2020) arXiv:2002.05191 [hep-ph]. [48] B. Schenke, C. Shen, and D. Teaney, (2020),arXiv:2004.00690 [nucl-th].[49] S. R. de Groot, W. A. van Leewen, and C. G. van Weert, Relativistic Kinetic Theory: principles and applications (Elsevier North-Holland, 1980).[50] H. Grad, Commun.Pure Appl.Math. , 331 (1949).[51] A. H. Mueller, Phys.Lett. B475 , 220 (2000), arXiv:hep-ph/9909388 [hep-ph].[52] S. de Groot, W. van Leeuwen, and C. van Weert,
Relativistic Kinetic Theory: Principles and Applications (North Holland, Amsterdam, 1980).[53] A. Behtash, C. N. Cruz-Camacho, S. Kamata, andM. Martinez, Phys. Lett.
B797 , 134914 (2019),arXiv:1805.07881 [hep-th].[54] A. Behtash, S. Kamata, M. Martinez, and H. Shi, Phys.Rev.
D99 , 116012 (2019), arXiv:1901.08632 [hep-th].[55] A. Behtash, S. Kamata, M. Martinez, and H. Shi,(2019), arXiv:1911.06406 [hep-th].[56] M. Strickland, JHEP , 128 (2018), arXiv:1809.01200[nucl-th]. [57] A. Dash and V. Roy, (2020), arXiv:2001.10756 [nucl-th].[58] A. Kurkela and A. Mazeliauskas, Phys. Rev. Lett. ,142301 (2019), arXiv:1811.03040 [hep-ph].[59] K. Dusling, G. D. Moore, and D. Teaney, Phys. Rev. C81 , 034907 (2010), arXiv:0909.0754 [nucl-th].[60] R. Baier, P. Romatschke, D. T. Son, A. O. Starinets,and M. A. Stephanov, JHEP , 100 (2008),arXiv:0712.2451 [hep-th].[61] D. Teaney and L. Yan, Phys. Rev.
C89 , 014901 (2014),arXiv:1304.3753 [nucl-th].[62] A. Jaiswal, Phys. Rev.
C87 , 051901 (2013),arXiv:1302.6311 [nucl-th].[63] A. Kurkela, A. Mazeliauskas, J.-F. Paquet, S. Schlicht-ing, and D. Teaney, “KMPST: linearized kinetic theorypropagator of initial conditions for heavy ion collisions,”(2018), https://github.com/KMPST/KoMPoST.[64] G. S. Denicol and J. Noronha, in28th InternationalConference on Ultrarelativistic Nucleus-Nucleus Colli-sions (Quark Matter 2019) Wuhan, China, November 4-9, 2019