Hydrodynamization and non-equilibrium Green's functions in kinetic theory

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 and to what extent a hydrodynamic description of the system becomes applicable. Within the kinetic theory framework we derive a new method to calculate time dependent non-equilibrium Green's functions describing the evolution of energy and momentum perturbations on top of an evolving far-from-equilibrium background. We discuss the approach towards viscous hydrodynamics along with the emergence of various scaling phenomena for conformal systems. By comparing our results obtained 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 are sensitive to the underlying microscopic dynamics.


I. INTRODUCTION
Ultrarelativistic heavy ion experiments carried out at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) have produced a new state of matter where quarks and gluons are liberated from the incoming nuclei [1][2][3][4] Since the lifetime of this Quark-Gluon Plasma (QGP) is very short, its properties are reconstructed by analyzing a large set of hadronic observables. The phenomenological studies of the collected wealth of experimental data have shown that hydrodynamical models provide robust tools to explain the dynamics of the QGP [5][6][7][8][9]. As a result, a new paradigm in high energy nuclear physics has emerged where hydrodynamics plays a central role.
One of the most consistent findings in the hydrodynamical modeling of heavy ion collisions is a small value of the shear viscosity over entropy ratio η/s ∼ 1/(4π), extracted from large number of observables measured in ultrarelativistic heavy ion collisions over a different range of energies at both RHIC and LHC. Currently, the standard approach to describe high-energy heavy ion collisions is based on multistage models where the space-time dynamics of the QGP is described using relativistic viscous hydrodynamics. Subsequently, the energy density of the QGP is converted into a hadron gas which eventually freezes out, yielding the final state hadronic observables that can compared to experimental measurements. One key ingredient in this phenomenological approach are the initial conditions, characterizing the initial distributions of the energy density and flow velocities, which are then * skamata11phys@gmail.com † mmarti11@ncsu.edu ‡ p.plaschke@uni-bielefeld.de § s.ochsenfeld@uni-bielefeld.de ¶ sschlichting@physik.uni-bielefeld.de propagated via hydrodynamics. Clearly, the calculation of these quantities from first principles QCD represents an enormous challenge, which despite important developments during the last years [10][11][12][13][14][15][16][17][18] has not been answered completely. One important feature shared by all initial state models is related to the far-from-equilibrium nature of the QCD matter, produced immediately after the collision of heavy nuclei on a time scale τ 0 1fm/c [19]. Since at these early times the QCD matter is also subject to a rapid longitudinal expansion, viscous hydrodynamics which is an effective theory for the long-time and long wave-length behavior close to equilibrium is not necessarily applicable at such early times. Although in recent years the existence of a far-from-equilibrium fluid dynamical theory has been advocated , its formulation remains to be completed and thus, in practice hydrodynamic simulations of heavy-ion collisions are usually initialized after a certain period time τ = τ hydro ∼ 1fm/c, where the long. expansion is less rapid and the QGP has evolved towards local thermal equilibrium.
The question how to match the non-equilibratium initial state at τ 0 to the initial state for hydrodynamics at τ hydro was the subject of a series of recent papers [43][44][45][46]. Based on an underlying microscopic description in QCD kinetic theory, the initial conditions for the energy-momentum tensor T µν (τ hydro ) at the beginning of the hydrodynamic phase are hereby obtained from the initial energy-momentum tensor T µν (τ 0 ), by determining the evolution of macroscopic quantities, based on non-equilibrium Green's functions of the energy momentum tensor [44,45]. This new framework, dubbed as KøMPøST [44,45], has proven to be a powerful tool to describe the pre-equilibrium evolution of heavy-ion collisions on an event-by-event basis [47,48]. While the detailed phenomenological consequences of the preequilibrium stage remain to be explored, it has been demonstrated that the subsequent hydrodynamic evolu-tion becomes independent of the hydrodynamic initialization time τ hydro as long as the latter is chosen to be in the regime where both kinetic theory and hydrodynamic descriptions overlap.
While [44,45] presented a calculation of nonequilibrium Green's functions in QCD kinetic theory, it is also interesting to explore to what extent the microscopic details affect the evolution of the macroscopic quantities far-from-equilibrium. In this work we therefore present a new method to calculate non-equilibrium Green's functions of the energy momentum tensor. We generalize the method of moments approach [49,50] by incorporating the response of a far-from-equilibrium expanding plasma against linear perturbations, and analyze their non-equilibrium in the relaxation time approximation. Instead of solving linear kinetic theory for perturbations of the phase-space distribution function we analyze the equations of motion of the corresponding linearized moments. Eventually, the Green's functions of the energy momentum tensor are reconstructed from the linearized moments.
This paper is organized as follows: In Sect. II we introduce the relevant aspects of the Boltzmann equation. We explain the main aspects of the emergent attracting behavior of the non-equilibrated Bjorken flow background in Sect. III. In Sect. IV we present the formalism to describe the space-time evolution of the perturbations. Using this formalism we proceed in Sect. V to calculate the Green's functions of the energy-momentum tensor. Conclusions and outlook are discussed in Sect. VI. Some of technical aspects of our work are briefly discussed in Appendix A.

II. BOLTZMANN EQUATION WITHIN THE RELAXATION TIME APPROXIMATION
Starting point of our analysis is the Boltzmann equation within the relaxation time approximation (RTA) where the coordinate system defined in Minkowski with the local rest-frame velocity u µ (x) determined via the Landau matching condition The fluid velocity is defined as a time-like eigenvector (u 2 = +1) of the energy momentum tensor for on-shell particles where we denote the on-shell momentum average of any observable as The effective temperature T (x) entering in Eq. (2b) is determined from the (equilibrium) equation of state e eq (T ), by matching the corresponding eigenvalue e to the equilibrium energy density. If not stated otherwise, we will consider an ultra-relativistic system of massless bosons where the equilibrium distribution function is given by the Bose-Einstein distribution f eq. (x) = 1/ (e x − 1), such that We are interested in longitudinal boost-invariant expanding system and thus, we use the hyperbolic Bjorken coordinates defined in terms of the cartesian coordinates as so the metric g µν = diag(1, −1, −1, −τ 2 ) and − g(x) = τ . Similarly, the four-momentum of a relativistic massless particle is with y = arctanh (p 3 /p 0 ) and p T = |p|.
In the Bjorken coordinates the Boltzmann equation (1) takes the following form where Hereafter we shall use the roman letter i = x, y to denote the summation only over the transverse coordinates. By virtue of the coordinate transformation (5), the derivatives w.r.t. τ, ς in Eq. (7) are taken at constant p T and momentum space rapidity y. However, when analyzing the dynamics of a boost invariant medium it is more convenient to work with the (dimensionless) longitudinal momentum variable By transforming the space-time derivatives according to one finds that the two additional terms cancel each other, such that which is the form of the equation that we will consider in this work. We note in passing, that it is also common in the literature to express the dynamics in terms of the longitudinal momentum in the local rest frame In this case only the derivative w.r.t. the longitudinal rapidity is affected by the transformation and the Boltzmann equation takes the form (see e.g. [51])

III. EVOLUTION OF BOOST INVARIANT HOMOGENOUS BACKGROUND
During the early pre-equilibrium stage of a heavy-ion collision, which lasts about 1fm/c, the non-equilibrium plasma is subject to a rapid longitudinal expansion. Conversely, in the transverse plane the plasma is initially created at rest, and the transverse expansion only builds up in response to local energy-density gradients on a time scale τ ∼ R, where R denotes the system size. Due to this separation of scales, it is a reasonable assumption to neglect the transverse expansion during the pre-equilibrium stage, and first consider an idealized situation of Bjorken flow, where the system is longitudinally boost-invariant, parity invariant under spatial reflexions along the longitudinal beam line and azimuthally symmetric and translationally invariant in the transverse plane. Space-time dependent variations can subsequently be addressed by studying small deviations from this average background behavior, and will be discussed in Sec. IV.
In this section we discuss the space-time evolution of an expanding background undergoing Bjorken expansion, such that the aforementioned symmetries constrain the functional form of the distribution function for the background to be of the form The energy momentum tensor has only non-vanishing diagonal components, i.e., with the energy density (e), transverse and longitudinal pressures (p T /L ) determined by Due to scale invariance, the previous expressions automatically satisfy the tracelessness condition e = 2 p T +p L . Based on the explicit form of T µν BG the Landau matching condition in Eq. (2) becomes trivial with and the kinetic equation for the evolution of the background distribution takes the familiar form where p τ = p 2 T + (p ς /τ ) 2 according to the on-shell mass condition.

A. Evolution equations for moments
Several strategies have been explored in the literature to solve Eq. (20) [52]. Here we follow Grad's original approach [50] where instead of finding solutions for the phase-space distribution function f (x, p) we study the dynamics of its moments. The latter are directly connected to macroscopic observables as we indicate below.
We consider the following moments C m l (τ ) = ν eff dp ς (2π) where the angles are defined as cos θ p = p ς /(τ p τ ), tan φ p = p 1 /p 2 in the co-moving coordinates and Y m l denote the spherical harmonics with normalization Non-vanishing components of the background energymomentum tensor in (17) are related with the moments C m l in (21) as follows Specifcially, if the background distribution function is in thermal equilibrium, one has C m l | eq (τ ) = where the Landau matching conditions (2) was enforced. Now the evolution equations of those moments is simply obtained by taking the explicit time derivative in its definition. After some careful algebra and using a series of identities of the spherical harmonics (see App. A), the evolution equation for the moments C m l takes the following form where the coefficients b m l,−2 , b m l,0 and b m l,+2 are given by We note that in Eq. (27) only the moments l and l ± 2 are coupled reflecting the parity symmetry of the background. Now due to the azimuthal symmetry of the background, the evolution of moments with different m are not coupled and only moments with m = 0 are nonvanishing. Furthermore, Eqs. (27) present interesting mathematical resurgent properties which were discussed extensively in Refs. [53][54][55].

B. Initial conditions
Since at early times τ τ R the system is unable to sustain sizeable longitudinal momenta, the physically relevant initial conditions for the phase-space distribution are naturally of a form where p p T . Indeed, previous works [35-37, 53, 55-57] have shown that the extreme limit where the (longitudinal) support of the phase-space distribution shrinks to a Dirac delta function corresponds to a non-equilibrium attractor of the kinetic equation. We will therefore consider precisely this case and choose the initial phase-space distribution as where the normalization chosen such that the initial energy density per unit rapidity remains constant, i.e.
By construction, the initial condition in Eq. (29) fixes the initial values of the moments C m l at τ 0 , i.e., such that for m = 0 one has Interestingly, one finds that at the level of relaxation time approximation, the shape of the azimuthally symmetric momentum distribution dN dςd 2 pd 2 x is completely irrelevant for the dynamics. Noteably, this is in sharp contrast to more realistic description in Yang-Mills kinetic theory [23,44,45,58], where different processes (e.g. inelastic and elastic interactions) exhibit different parametric dependencies on the momenta (see e.g. [19,59]).

C. Evolution of energy momentum tensor for constant relaxation time (τR = const)
We first analyze the evolution of the energy momentum tensor for a constant relaxation time τ R = const. By truncating the evolution equations up to a certain value of l < l max , setting C m l = C m l | eq = 0, we can obtain a numerical solution of the coupled set of evolution equation as a function of the dimensionless evolution time variable τ /τ R . Noteably the truncation scheme converges rapidly except for very small values of τ /τ R 1 where higher moments play an important role. However, for the typical regime of interest no visible deviations of the results for the lowest l ≤ 2 moments occur for l max 32, and if not stated otherwise we use l max = 512 to produce the figures.

Comparison with hydrodynamics
Our results are compactly summarized in Fig. 1, which shows the non-equilibrium attractor solution for the evolution of the energy and pressure densities. For the Bjorken flow we can also compare the numerical solution of the Boltzmann equation in relaxation time approximation to the truncation in relativistic viscous hydrodynamics. Starting from the conservation equation the longitudinal and transverse pressure within the hydrodynamic approach are given by 3p T e 3p L Energy/pressure density: Navier-Stokes Hydro where π ς ς is the only independent component of the shear viscous tensor. Now π ς ς is expanded up to second order in the gradient hydrodynamical expansion as [60] where for the system under consideration the equation of state and transport coefficients in the RTA approximation are determined by [37,54,61,62] p = e/3 , such that η/τ = 4 15 τ R τ e. Hence the asymptotic solution for the energy density up to second order within the gradient hydrodynamic expansion (or alternatively when τ /τ R 1) takes the form with the asymptotic integration constant τ 4/3 e ∞ determined as from the numerical solution of the Boltzmann equation in relaxation time approximation. By comparing the different curves in Fig. 1 one observes that viscous hydrodynamics starts to describe the evolution of the energy momentum tensor around evolution times τ /τ R ∼ 2 − 3. Since including the second order correction does not significantly improve the agreement in Fig. 1, we only present the Navier-Stokes limit, i.e. the leading term in the asymptotic expansion in τ R /τ in Eq. (38).
Noteably, one can also define an effective ratio between the inverse Reynolds and the Knudsen number as follows with Kn = τ R /τ being the Knudsen number and the inverse Reynolds number Re −1 is Near equilibrium where the gradient expansion (36) holds, one can express the inverse Reynolds number as While the gradient expansion (36) breaks down at early times, where the Knudsen number Kn 1 , the ratio in Eq. (40) can still be used to quantify the evolution. Specficially, for the case of a constant relaxation time, one can use Eq. (37) to define an effective relaxation rate τ ef f R away from equilibrium which is constructed such that at late times the ratio in Eq. (43) converges to unity as expected.
Numerical results for Fig. 1 indicate that at early times where the Knudsen number Kn 1, the effective relaxation rate is significantly reduced by the inclusion of higher order dynamical moments C m l . By explicitly comparing different truncations (l max = 4, 8, · · · ) of the infinite hierarchy of moment equations, one also observes a rapid convergence in the sense that the evolution of the low order moments becomes increasingly insensitive to the higher order moments.  We now investigate the more commonly studied case where the relaxation time is chosen inversely proportional to 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 timelike variable which one can also identify as the inverse Knudsen number x ≡ Kn −1 . Since for the conformal system, the relaxation time depends on the temperature of the system as τ R = 5η/sT (τ ), it is convenient to perform a change of variables to express where we defined the scale-factor and used the equation of motion for the energy density along with the fact that for an ultra-relativistic system dT T = de 4e . The evolution equation for the moments can be then re-cast into the form indicating that for a given initial condition the nonequilibrium evolution of the system is uniquely determined by the conformal scaling variable x.
Numerical results for the evolution of the background energy-momentum tensor are presented in Fig. 2, where we also compare to the results for Yang-Mills kinetic theory obtained in [44,45]. Clearly the overall behavior of the curves is quite similar, showing a smooth transition from an approximate free-streaming behavior at early times towards the universal hydrodynamic behavior at late times. Similarly, to our previous discussion, this behavior can also be analyzed in terms of an effective shear viscosity over entropy ratio in far from equilibrium regime, which for the conformal system is given by which is again constructed such that (η/s) ef f η/s approaches unity in the limit τ /τ R 1. Even though this ratio depicted in the lower panel of Fig. 2 magnifies the differences between the results for Yang-Mills kinetic theory [44,45] and the relaxation time approximation, the overall differences are at most at the 10% level for intermediate timesw ∼ 1.
Nevertheless, the fact that the approach towards visc. hydrodynamics is different for Yang-Mills kinetic theory and RTA results in a mismatch in the ratios of the initial energy density to final energy density. Since the early time behavior is governed by free-streaming, one finds that e(τ ) ∝ 1/τ for x 1 such that the attractor curve can be parametrized as By inverting this relation, one then obtains the energy density at late times as [46] Specifically, for the conformal relaxation time approximation, we find C ∞ ≈ 0.9 whereas for the Yang-Mills kinetic theory results of [44,45] the pre-factor C ∞ ≈ 1 is about ten percent larger [46].

IV. ENERGY MOMENTUM PERTURBATIONS AROUND BJORKEN FLOW
So far we have addressed the non-equilibrium evolution of the average boost invariant and homogenous background. We will now consider the propagation of lin-earized perturbations, sourced by (small) space-time dependent deviations of the initial energy momentum tensor from its (local) average. By linearizing the kinetic equation around the boost invariant and homogenous background one finds an evolution equation for the perturbation of the distribution function δf , i.e., where f eq (x) = df eq (x)/dx denotes the derivative of the equilibrium distribution. In the above expression, the change in the rest-frame velocity δu µ (x) and local equilibrium temperature δT (x) are to be determined selfconsistently from the (linearized) Landau matching condition. Starting from the linearized perturbations of the energy momentum tensor the change in the rest-frame velocity δu µ and energy density in the local rest-frame δe are determined from the linearized eigenvalue equation with u µ δu µ = 0 .
By using the leading order solution u µ T µν = eu ν of the eigenvalue problem, one reads off Our strategy to determine the evolution of energy and momentum perturbations then consist in determining all the coefficients on the r.h.s. of Eq. (52) based on certain moments of the perturbation of the phase-space distribution function. In the following section we shall explain this procedure in detail.
A. Evolution equations for energy-momentum perturbations in the transverse plane We will from now on restrict our attention to energymomentum pertubations in the transverse plane, i.e. we will only consider variations in the transverse coordinates x. Since we are considering the evolution of linearized perturbations on top of a (transversely) homogenous background, it is natural to express them in a Fourier basis such that for each k-mode where we denote δf k (τ, p, |p ς |) ≡ δf (τ, k, p, |p ς |). In our approach the Z 2 subgroup symmetry of reflections along the beam axis is not broken at the level of the perturbations so the longitudinal rest frame velocity δu ς vanishes identically. The transverse flow velocity can be decomposed in the components parallel (δu k (τ )) and perpendicular (δu ⊥ k (τ )) to the wave vector k. Hence the perturbations to the thermodynamic fields take the following form where for physical perturbations, the reality conditions of the perturbations imply δT k (τ ) = δT * −k (τ ) and δu i with φ pk = φ p − φ k being the angle between p and k in the transverse plane and sin θ p = |p|/p τ and inserting the above expressions into the linearized kinetic equation (52) one then finds where the terms on the left hand side describe freestreaming, the terms in the first line correspond to the relaxation of the perturbation, the terms in the second line describe the change of the equilibrium distribution due to the perturbations and the terms in the last line describe the change in the relaxation of the non-equilibrium background due to perturbations.

B. Evolution equation of the spherical harmonic moments
We follow the same strategy as for the evolution of the background and transform the evolution equation for the distribution function into a coupled set of evolution equations for the spherical harmonic moments δC m l,k (τ ) = ν eff dp ς (2π) where the azimuthal and polar angles φ pk , θ p are measured with respect to the transverse wave vector (k) and the longitudinal rapidity (ς) axis. 1 Similarly as for the background, the various components of the energymomentum tensor are explicitly given in terms of the lowest order ( = 0, 1, 2) moments, as where we have conveniently decomposed all transverse vectors δT τ i , δT ςi into components parallel, perpendicular and independent w.r.t. to the wave-vector k/|k| (and similarly for the tensor δT ij ). Due to the residual Z 2 symmetry of long. reflection in rapidity, the components T τ ς and T ςi vanish identically, whereas all other components can in principle be non-zero in our setup. Evaluating the various couplings between spherical harmonics using the identities listed in App. A, the evolution equations for the spherical harmonic moments δC m l then take the form 1 Since physical perturbations of the phase-space distributions δf (τ, x, p, |pς |) are real-valued, the linearized perturbations in Fourier should also satisfy the condition δf −k (τ, p, |pς |) = δf * k (τ, p, |pς |), such that in terms of the spherical harmonic moments one finds δC m ,−k (τ ) = (−1) m δC −m ,k (τ ) * for physical perturbations.
τ ∂ τ δC m l,k = b m l,−2 δC m l−2,k + b m l,0 δC m l,k + b m l,+2 δC m l+2,k − i|k|τ 2 u m l,− δC m+1 l−1,k + u m l,+ δC m+1 l+1,k + d m l,− δC m−1 l−1,k + d m l,+ δC m−1 where the coefficients b m l are the same as in Eq. (28) and Physically the terms in the first line of Eq. (64) correspond to a free-streaming evolution, while the terms in the last few lines capture the relaxation towards equilibrium, including the changes of the equilibrium distribution and background equilibration. By C eq in Eq. (64) we denote the moments (C eq ) m l (τ ) = dp ς (2π) which can be determined via integration by parts as with the equilibrium moments (C eq ) m l determined by Eq. (26). By inserting the relations in Eqns. (63) into the linearized Landau matching conditions in Eq. (56), one also finds that the linearized perturbations of the energy density δe k and long. and transverse flow velocities δu k and δu ⊥ k are given by such that the evolution equations for the moments in Eq. (64) form a closed set of equations. We further note that by virtue of the decomposition into spherical harmonic moments, Y m l (φ pk , θ p ), the information on the direction of the transverse wave-vector k has disappeared from the evolution equation, which as a consequence of the azimuthal rotation symmetry of the background only depends on the magnitude of the wave-vector |k|.
Before we address the physically relevant initial conditions for the moments δC m l,k , we also note that it is often useful to consider the evolution equation at a fixed value of propagation phase rather than a fixed value for the wave-number |k|. By changing the variables from |k| to κ = |k|(τ − τ 0 ) for the moments, the time derivate needs to be evaluated according to such that the evolution equation for the moments receives one additional term associated with this change of variables. Similarly, for a conformal system, it is also convenient to express the evolution in terms of the scaled evolution time x = τ /τ R introduced in Eq. (44). Starting from Eq. (64) and denoting s(τ ) = τ −τ0 τ , with to perform these changes, the equation of motion for the moments then takes the form which we will employ below to obtain the numerical solution for the evolution of perturbations.
Based on Eq. (73), one also explicitly observes that in the limit τ 0 /τ R → 0 where s(x) = 1 is a fixed point of Eq. (72), the solution to the evolution equation for the moments δC m l,k only depends on the propagation phase κ = |k|(τ − τ 0 ) and the scaled evolution time x = τ /τ R = T (τ )τ 5η/s . While this conformal scaling behavior was empirically observed in [44,45] from numerical solutions of the Boltzmann equation in QCD kinetic theory at different values of the coupling strength λ = g 2 N c , it is interesting to point out that in the present context the conformal scaling behavior directly manifests itself at the level of the equations of motion.

C. Initial conditions for energy and momentum perturbations
So far we have discussed, the evolution equations for linearized perturbations on top of a Bjorken background. Now in order to apply this framework to describe the early time dynamics of high-energy heavy-ion collisions, the equations of motion need to be supplemented by suitable initial conditions, which describe the associated change of the phase-space density at initial time. While in principle one could imagine a large variety of different initial conditions, we will follow [44,45] and only consider the response of the system to changes of the conserved quantities of the system, associated with initial energy and momentum perturbations as detailed below.

Energy perturbations
We follow the arguments of [44] and associate initial energy perturbations with an infinitesimal change of the energy scale of the background distribution, such that the associated phase-space distribution for energy per-turbations is given by where we introduced the phase-factor e −ik p |p| τ0 to account for the free-streaming evolution at early times τ < τ 0 τ R . Based on the explicit form of the initial background distribution in Eq. (29), the integrals for the moments δC m l in Eq. (62) can be evaluated using 1 2π where (eτ ) 0 denotes the asymptotic energy density of the background (c.f. Sec. III). Specifically, for the energy and momentum density one has reproducing the result of [44] for the free-streaming response function.

Momentum perturbations
Similarly, we associate initial momentum perturbations with an infinitesimal change of the transverse velocity of the background, such that following the arguments of [44] the associated phase-space distribution for momentum perturbations is given by where the index i contains the information about the direction of the initial momentum perturbation. Decomposing δf i k into the directions parallel and perpendicular to the wave-vector, we can distinguish between longitudinal and transverse momentum perturbations Evaluating the moments, one then finds for longitudinal momentum perturbations and similarly for transverse momentum perturbations Specifically, for the energy and momentum density one has for longitudinal momentum perturbations, whereas for transverse momentum perturbations in agreement with [44].

V. NON-EQUILIBRIUM GREEN'S FUNCTIONS OF THE ENERGY-MOMENTUM TENSOR
We now proceed to the calculation of the response of the energy-moment tensor to initial energy and momentum perturbations, based on numerical solutions of the evolution equations for the moments δC m l,k starting from the initial condition described in the previous section. Since the set of moments δC m l,k contain an overwhelming amount of information, we will therefore restrict our attention to the evolution of the low order moments δC m l,k with ≤ 2, which can be related to the various components of the energy momentum tensor δT µν k according to Eq. (63). Instead of investigating the dynamics of individual moments δC m l,k directly, we find it more insightful to consider the linear response functions G µν αβ (k, τ, τ 0 ) introduced in [44,45], such that Noteably, the Green's functions G µν αβ (k, τ, τ 0 ) provide the builiding block of the pre-equilibrium computer code KøMPøST [63] and can be obtained in terms of linear combinations of the moments as described below. Since we are primarily interested in the limit τ 0 /τ R → 0, where the kinetic theory framework describes the equilibration process of the system all the way from very early times up to the onset of hydrodynamic behavior, we will drop the explicit dependence on τ 0 in the following to lighten the notation. By expressing the response to an initial energy perturbation in the following basis of scalars (s), vectors (v) and tensors (t) and adapting the normalization δe(τ 0 )/e(τ 0 ) = 1, the relevant response functions can then be determined from (c.f. Appendix B of [44]) Similarly, the response to an initial momentum perturbation can be characterized by a set of six independent response functions,

A. Numerical results
We will now present numerical results for the various response functions, focusing on the case of a conformal relaxation time. We follow the same strategy as for the background and solve a truncated set of evolution equa-tions, checking that including higher order moments does not significantly alter the results. Noteably, we find that a rather large set of moments is required to achieve apparent convergence and we will present results for l max = 64 in the following.
Our results for the evolution of the response func- tions are compactly summarized in Figs. 3 and 4 , where we present the spectrum of perturbations G µν τ τ (k, τ − τ 0 , τ /τ R ) and G µν τ i (k, τ − τ 0 , τ /τ R ) as a function of |k|(τ − τ 0 ) for different values of the evolution timẽ w = T (τ )τ /(4πη/s). Different panels in Figs. 3 and 4 show the results for the different (µν) components, decomposed in the tensor basis described above. In order to facilitate the interpretation, we have also labeled the various response functions according to the components of the energy-momentum tensor that they affect.
Starting from the initial free-streaming behavior at early times T (τ )τ /(4πη/s) 1 discussed in Sec. IV C, one observes how towards later times the viscous damping of large |k|(τ − τ 0 ) modes sets in, such that by the time the system enters the hydrodynamic regime (w = T (τ )τ /(4πη/s) ∼ 1) only the long wave-length modes survive. Since at early times the system is highly anisotropic, the longitudinal pressure is effectively zero and transverse perturbations initially propagate with a phase-velocity of nearly the speed of light. Subsequently, as the system becomes more and more isotropic the phase-velocity decreases and eventually approaching the speed of sound, which in Figs. 3 and 4 results in a shift of the peaks towards larger values of the propagation phase |k|(τ − τ 0 ). Strikingly, the qualitative behavior observed from Figs. 3 and 4 is very similar to the results obtained in Yang-Mills kinetic theory in [44], albeit we find that in the relaxation time approximation the viscous damping of short wave length modes becomes efficient on a somewhat shorter time scale.
Even though some of the features of the evolution can be understood quite naturally in wave-number (k) space, in practice one is mostly interested in the Green's functions G µν τ τ (x − x 0 , τ − τ 0 , τ /τ R ) in position space, which directly describe the physical response of the energy momentum tensor δT µν (τ, x) to a localized initial energy perturbation δT τ τ (τ 0 , x 0 ) according to  6. Evolution of the energy-momentum response to an initial momentum 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 function of the propagation distance |x − x0|/(τ − τ0).
In practice, the coordinate space response can be obtained in a straightforward way via a set of Bessel-Fourier transforms of the response functions, and we use 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 results in [44] are presented in terms of the time variable x Id S = T Id (τ )τ /(4πη/s), where T Id (τ ) is defined via the asymptotic temperature according to T Id (τ ) = τ −1/3 lim τ →∞ (T τ 1/3 ). However, from the point of view of the RTA Boltzmann equation it is more natural to study the evolution as a function of the scaling variablẽ w = T (τ )τ /(4πη/s), where T (τ ) denotes the equilibrium temperature obtained from e(τ ) via Landau matching. Of course, the two quantities are related by x Id S = w (τ 4/3 e)∞ τ 4/3 e(τ ) 1/4 and we have taken this difference into account in all explicit comparisons presented in this paper. In order to provide an apples to apples comparison between results obtained within the relaxation time approximation and Yang-Mills kinetic theory, the different response functions are smeared out with the same smearing kernel K σ = exp(−k 2 (τ − τ 0 ) 2 /2σ 2 ) as in [44], which for Yang-Mills kinetic theory results is necessary in order to stabilize the numerical Bessel-Fourier transform.
Our results for the coordinate space response functions are shown in Figs. 5 and 6, where we display the various response functions as a function of the propagation distance |x−x 0 |/(τ −τ 0 ). Solid curves in Figs. 5 and 6 show our results obtained from the Boltzmann equation in relaxation time approximation, which are compared to the results obtained in Yang-Mills kinetic theory from [44] shown as dashed curves. Some of the most important features that can be immediately observed from Figs. 5 and 6, include the viscous broadening of the peaks due to the damping of high wave-number modes, as well as the shift of the peaks towards smaller values of |x − x 0 |/(τ − τ 0 ) associated with the aforementioned change of the effective (transverse) speed of sound, due to the increase of the longitudinal pressure. One also observes a decrease of the shear-stress response compared to the pressure response (for both energy and momentum perturbations), such that byw = T (τ )τ /(4πη/s) ∼ 1 when the background evolution starts to be captured by viscous hydrodynamics, also the dissipative corrections to the perturbations become sub-leading.
While the qualitative behavior of the various response functions is quite similar for RTA and Yang-Mills kinetic theory, one clearly observes that the RTA shows a faster departure from the early free streaming behavior, leading to differences ∼ 10% on the relevant time scales w = T (τ )τ /(4πη/s) ∼ 1 where viscous hydrodynamics becomes applicable. Even though these differences appear to be rather small, it would nevertheless be interesting to explore to what extent such differences in the early-time non-equilibrium dynamics can manifest themselves in final state observables in high-energy heavy-ion collisions.

VI. CONCLUSIONS & OUTLOOK
We derived a new method to calculate non-equilibrium Green's function of the energy momentum tensor based on moment equations of kinetic equations for linearized perturbations. Due to the particularly simple structure of the relaxation time approximation considered in this work, we obtained a closed set of moment equations for the evolution of the dimension four moments C m l and δC m l , which are relevant to study the evolution of the energy-momentum tensor, that is determined by the lowest order ( ≤ 2) moments. Even though, for more complex interactions, a truncation at the level of dimension four operators is no longer sufficient to obtain a closed set of evolution equations, we naturally believe that by including higher order operators our method can be extended to systems with more complex interactions, thereby generalizing the usual moment method at level of perturbations.
Based on the evolution equations for the moments C m l and δC m l in Eqns. (27) and (64), we studied the evolution of average energy-momentum tensor in Bjorken flow, as well as the out-of-equilibrium linear response of the system to initial energy and momentum perturbations in the transverse plane. By truncating the infinite hierarchy of moment equations at large finite order, we obtained numerical solutions for the evolution of the non-equilibrium background and the Green's functions. When comparing our results to previous calculations of KøMPøSTin Yang-Mills theory [44], we found a striking similarity between the different theories. Even though the macroscopic differences between the two microscopic calculations are only at the ten percent level, it would be interesting to explore to what extent these can affect concrete observables, such as e.g. the flow harmonics v n , the charged particle multiplicity dN ch /dη or the transverse energy dE ⊥ /dη. Since there are first hints that in particular the charged particle multiplicity dN ch /dη, may be a rather sensitive measure of the entropy production during the pre-equilibrium phase [46], this remains an interesting question which we expect to be addressed in more detail in future studies.
Besides our numerical studies, we also found that various conformal scaling features, which have been empirically observed in [44], can be directly seen at level of the equations of motion for the moments δC m l , and we expect that further analytic insights into the structure of the non-equilibrium Green's functions. Specifically, it would be interesting to further explore for example the early and late time asymptotics of the Green's functions based on this formulation. Furthermore, one can also investigate to what extent the highly non-trivial evolution of the Green's functions can be captured by "renormalized transport coefficients" as suggested in various works [36,53,64] studying the highly symmetric Bjorken flow.
Beyond such analytic insights, it would also interesting to further systematically extent the pre-equilibrium description in terms of non-equilibrium Green's functions. Since the numerical solution of the moment equations is comparatively straightforward, the Boltzmann equation in relaxation time approximation provides an ideal testing ground for such ideas. Specifically, with the methodology developed in this paper it should be comparatively straightforward to extent the description to include also longitudinal fluctuations, investigate Green's functions for higher order moments or study the effect of additional conserved charges, all of which are rather challenging tasks within a QCD kinetic description. By explicitly comparing the linearized description in terms of non-equilibrium Green's functions to full numerical solutions of the RTA Boltzmann equation, one could also obtain additional insights into the reliability and breakdown of the linearized description and potentially improve the range of applicablity to describe small collision systems. , which upon further simplifications yield the results quoted in the main text.