Effective kinetic description of event-by-event pre-equilibrium dynamics in high-energy heavy-ion collisions

We develop a macroscopic description of the space-time evolution of the energy-momentum tensor during the pre-equilibrium stage of a high-energy heavy-ion collision. Based on a weak coupling effective kinetic description of the microscopic equilibration process (\`a la"bottom-up"), we calculate the non-equilibrium evolution of the local background energy-momentum tensor as well as the non-equilibrium linear response to transverse energy and momentum perturbations for realistic boost-invariant initial conditions for heavy ion collisions. We demonstrate how this framework can be used on an event-by-event basis to propagate the energy momentum tensor from far-from-equilibrium initial state models, e.g. IP-Glasma, to the time $\tau_\text{hydro}$ when the system is well described by relativistic viscous hydrodynamics. The subsequent hydrodynamic evolution becomes essentially independent of the hydrodynamic initialization time $\tau_\text{hydro}$ as long as $\tau_\text{hydro}$ is chosen in an appropriate range where both kinetic and hydrodynamic descriptions overlap. We find that for $\sqrt{s_{NN}}=2.76\,\text{TeV}$ central Pb-Pb collisions, the typical time scale when viscous hydrodynamics with shear viscosity over entropy ratio $\eta/s=0.16$ becomes applicable is $\tau_\text{hydro}\sim 1\,\text{fm/c}$ after the collision.

High-energy heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) probe nuclear matter at extreme densities and temperatures. One of the primary goals of heavy-ion research is to study the properties of the new phase of deconfined matter created in such collisions: the Quark-Gluon Plasma (QGP).
Sophisticated multistage models, in which the QGP evolution is described by viscous relativistic hydrodynamics, have been remarkably successful in describing the soft hadronic observables measured in heavy-ion collisions [1][2][3][4][5]. Comprehensive model-to-data comparisons have been made to quantify systematically the constraints provided by measurements on the transport coefficient of the QGP, and to understand the impact of different observables on these constraints [6][7][8].
Considerable progress has been made to increase the predictive power of hydrodynamic simulations and to fully understand the assumptions built into such models. These advances include an evolving understanding of the conditions necessary for hydrodynamics to be applicable [9][10][11][12][13], and a more consistent treatment of the transition from hydrodynamics to hadronic kinetics at late times [14]. Progress is also being made on understanding the kinetics of early stages of the collision and the subsequent transition to hydrodynamics, which is the topic of this work.
The initial conditions of hydrodynamic models remain one of the major sources of uncertainty in phenomenological studies of heavy ion collisions. We will provide a practical way to propagate the energy-momentum tensor in a far-from-equilibrium initial state to a time when viscous hydrodynamics becomes applicable. Our goal is to have consistent overlapping descriptions of the early time dynamics, and to limit the dependence of the hydrodynamic model on ad hoc parameters such as the hydrodynamic initialization time τ hydro [15,16].
One approach to the initial conditions is simply to parameterize the initial energy density and its fluctuations, sidestepping the thermalization process with additional parameters. Glauber-based models are commonly used for this purpose, and provide the energy density at τ hydro [17][18][19][20]. Besides the energy density, the remaining hydrodynamic fields (such as the flow velocity and the shear and bulk tensors) must also be parameterized, leading to an uncomfortable growth in the number of free parameters. Physically motivated models such as free streaming [21,22] and the gradient expansion [15,23,24] have been used to relate the initial energy density to the full stress tensor which is ultimately needed to start the hydrodynamic simulation.
At weak coupling significant progress has been made in constructing a complete picture of the early time dynamics before τ hydro . Since the relevant degrees of freedom change as a function of time, a consistent theoretical description of the pre-equilibrium stages requires a combination of different weak-coupling methods. Based on the Color Glass Condensate framework (CGC), the initial state immediately after the collision is characterized by strong color fields, whose dynamics is essentially non-perturbative and best-described in terms of classicalstatistical field theory [25][26][27][28][29]. After a short period of time ∼ 1/Q s the system becomes increasingly dilute. Genuine quantum effects can then be no longer neglected, and the subsequent dynamics is better described in terms of QCD effective kinetic theory [30]. Several studies (which include some of the present authors) have investigated the various stages of the equilibration process in detail, including the early time dynamics using classicalstatistical real-time lattice techniques [31][32][33][34], as well as the subsequent approach towards local thermal equilibrium using effective kinetic theory simulations [35][36][37][38], i.e. the "bottom-up" thermalization scenario [39].
Although the output of classical field simulations has been used to initialize hydrodynamic codes [40], a consistent treatment at weak coupling would pass the classical output through the kinetic theory simulation to determine the initial conditions for the subsequent hydrodynamic evolution. In this paper we provide a concrete realization of this set of steps, allowing for an eventby-event description of the early time dynamics of high-energy heavy-ion collisions which smoothly approaches hydrodynamics.
Elaborating on the ideas formulated in Ref. [38], we describe the pre-equilibrium dynamics macroscopically in terms of non-equilibrium response functions of the energy-momentum tensor. Specifically, linearized energy and momentum perturbations are propagated on top of a boost-invariant and locally homogeneous background, and the energy-momentum tensor is evolved from a nonequilibrium initial state to a later time when viscous hydrodynamics becomes applicable. We demonstrate that the pre-equilibrium evolution smoothly matches onto hydrodynamics, and the subsequent hydrodynamic evolution becomes essentially independent of the matching time τ hydro .
In order to obtain a smooth transition from kinetic description to realistic viscous hydrodynamic evolution with typical shear viscosity over entropy ratio η/s ∼ 0.16, the coupling constant λ = N c g 2 (the single parameter of kinetic theory) has to be extrapolated to large values of λ = 10-25. For such values of λ, the entire nonequilibrium kinetic evolution is very well described by universal functions of scaled evolution time τ T /(η/s). The scalability of background and linear response functions greatly simplifies practical application, since the kinetic pre-equilibrium evolution needs only to be calculated once and then can be applied to any event-by-event hydrodynamic simulations. In order to facilitate the use of our results in phenomenological description of eventby-event high-energy heavy-ion collisions, we make public the linearized kinetic theory response functions and our implementation of the linear pre-equilibrium propagator KøMPøST [41].
The paper is organized as follows. We introduce a general macroscopic description of local pre-equilibrium evolution based on linear response theory out of equilibrium in Sec. II A, and discuss how to obtain the relevant inputs from an underlying microscopic description in effective kinetic theory in Sec. II B. We further study the equilibration of a locally uniform background in Sec. II C and derive local hydrodynamization time for realistic initial conditions in Sec. II D. In Sec. III A we provide the general decomposition of energy-momentum tensor response functions to initial energy and momentum perturbations and in Sec. III B we discuss their realization in effective kinetic theory. The implementation of the kinetic theory pre-equilibrium phase for hydrodynamic models of heavy-ion collisions is detailed in Sec. IV, and then applied to two types of initial conditions, MC-Glauber and IP-Glasma initial conditions, in Sections V A and V B respectively. We conclude with a compact summary of our findings and discussion of future directions in Sec. VII. Several appendices provide details on the background scaling functions (Appendix A), determination of the kinetic response functions (Appendix B), free streaming response functions (Appendix C), hydrodynamic response functions (Appendix D) and kinetic response in the low-k limit (Appendix E).

A. Macroscopic description of equilibration
Although the initial state shortly after the collision of two heavy nuclei is presumably very complicated, many of the microscopic details wash out over the first ∼ 1 fm/c as the Quark Gluon Plasma (QGP) equilibrates and becomes a hydrodynamically expanding fluid. Since the late time hydrodynamic behavior is fully characterized by the energy and momentum densities, it is conceivable that the most important features of the pre-equilibrium evolution also can be characterized by the energy-momentum tensor T µν . Following the computational strategy of Ref. [38], we will use QCD kinetics as a microscopic theory to determine the nonequilibrium evolution of T µν . This evolution propagates the non-equilibrium stress from an initial time τ EKT ∼ 0.1 fm to a time when hydrodynamics becomes applicable τ hydro ∼ 1 fm, and provides a map of the form T µν (τ EKT , x)| out-of-equilibrium −→ T µν (τ hydro , x), (1) as illustrated by Fig. 1.
We first note that, by virtue of causality, all contributions to the energy-momentum tensor at a given spacetime point (τ hydro , x) are fully determined by the initial conditions at earlier time τ EKT in the causal neighborhood of point x |x − x| < c(τ hydro − τ EKT ), (2) which is illustrated by a circle in Fig. 1. The energymomentum tensor can always be split into a locally homogeneous background and perturbations around it 1 T µν (τ EKT , x ) = T µν x (τ EKT ) + δT µν x (τ EKT , x ). (3) Since we anticipate the time scale of the equilibration process (τ hydro −τ EKT ) 1 fm to be small compared to the system size R Pb ∼ 5 fm, the long wavelength variations of the energy-momentum tensor within the causal circle are small ∼ (τ hydro − τ EKT )/R Pb 1. Short wavelength fluctuations (of order nucleon size or less) are 1/ √ N suppressed by the number of participant sources. Small perturbations around the background T µν x (τ EKT ) can be described (to first approximation) by linear response theory. Based on this approach, the energy-momentum tensor at (τ hydro , x) can be expressed as a sum of the evolved background T µν x (τ hydro ) and the response to the initial out-of-equilibrium perturbations δT αβ x : The first term represents the non-linear equilibration of the boost invariant homogeneous background. This contribution, discussed in detail in Sec. II C, has no transverse flow, but constitutes the major part of the final energy density and pressure. The second term in Eq. (4) is a convolution of initial perturbations with the response function G µν αβ (x, x 0 , τ hydro , τ EKT ), which in our work is the sole contributor to the off-diagonal components of T µν (τ hydro , x). The multicomponent structure of G µν αβ and its realization in kinetic theory is examined in Sec. III. Finally, note that the longitudinal expansion has been factored out from the response functions by normalizing perturbations with the background energy density.
By applying the propagation formula in Eq. (4) to all points x in the transverse plane, we construct initial conditions T µν (τ hydro , x), which can be used as input to the subsequent hydrodynamic description of high-energy heavy-ion collisions. Since the out-of-equilibrium evolution of the background and response functions smoothly approaches hydrodynamics, the initialization time τ hydro can be changed without affecting physical observables. Most importantly, this evolution consistently describes the pre-equilibrium production of entropy and transverse flow.
Eq. (4) provides a general framework to study the macroscopic features of the pre-equilibrium evolution. 1 We will consider a boost-invariant form of the energy-momentum tensor throughout this work and bold spatial vectors x lie entirely in the transverse plane. Note that in principle this framework could be extended to include an inhomogeneous background and variations of the energy-momentum tensor in rapidity.
However, its inputs, the non-equilibrium evolution of the background stress T µν x (τ ) and response functions G µν αβ (x, x , τ hydro , τ EKT ), have to be computed based on an underlying microscopic description. In this work we use a weakly coupled kinetic theory to compute these inputs, and extrapolate our results to realistic values of the coupling constant. However, we note that in future applications, these could be also based on different microscopic descriptions e.g. strongly coupled holographic models.
While the details of the calculation are described in the following sections, we point out an important feature of the dynamics that greatly simplifies the practical use of Eq. (4) for realistic values of coupling constant. Specifically, we find that for the relevant range of couplings λ (or equivalently η/s), all of the η/s dependence of the background and response functions is described by η/sindependent functions of a scaled time, τ T /(η/s). Thus, these universal functions can be tabulated for one value of η/s, and the results can be used in Eq. (4) to initialize hydrodynamic simulations over a significant η/s range. Interestingly, this scaling property also allows us to smoothly extrapolate the weakly coupled results into the strongly coupled regime [42,43].

B. Effective kinetic theory & bottom-up thermalization
In this work the microscopic description of equilibration is provided by QCD effective kinetic theory [30,37,38]. Specifically, for the boost-invariant SU(N c ) gluonic plasma considered here, the Boltzmann equation for the color and spin averaged gluon distribution function f x,p takes the form 2 C 2↔2 [f x,p ] is the collision integral for elastic scatterings at leading order in the coupling constant, λ = 4πα s N c . The elastic scattering matrix element, diverges at small momentum transfer, and is regulated by a screening mass which is adjusted so that the simulation reproduces the leading order drag and momentum diffusion coefficients for isotropic distributions [38,44]. Similarly, C 1↔2 [f x,p ] describes the leading order inelastic (particle number changing) bremsstrahlung processes. It includes the Landau-Pomeranchuk-Migdal suppression of FIG. 2. Different phases of weak coupling equilibration (à la "bottom-up") for a Bjorken expanding plasma. Different panels show the gluon momentum distribution functionf (τ, p) in the (p x , p z ) plane as a function of scaled time τ T id. /(4πη/s): (a) The initial distribution is broadened by elastic 2 ↔ 2 scattering; (b) inelastic 1 ↔ 2 splitting and mini-jet quenching build up a low-momentum thermal bath; (c) the distribution approaches local thermal equilibrium, but has significant corrections which are described by viscous hydrodynamics. The initial conditions forf are given by Eq. (7) and the evolution is done for the coupling constant λ = 10 (with η/s ≈ 0.62). The color scale represents the distribution function in the range 0.1 <f < 5.
collinear radiation, and also uses the isotropic screening approximation [37]. Details of the numerical implementation can be found along with detailed expressions for the matrix elements in Appendix A of Ref. [38].
We determine the microscopic input to Eq. (4) by performing two independent sets of calculations: one for the average background T µν described in Sec. II C, and one for the response functions G µν αβ described in Sec. III B. We calculate the evolution of the average energy-momentum tensor T µν by studying the equilibration of a spatially homogeneous "background" distributionf (τ, p), starting from an initial condition where we denote cos(θ) = p z /|p| and choose Q 0 = 1.8 Q s , ξ = 10, and A = 5.24 as in previous publications [37,38]. This specific form is motivated by classical simulations of the early time dynamics [45], where the phase space distribution at τ ∼ τ EKT approaches a scaling form characterized by a large initial momentum anisotropy and a transverse momentum scale ∼ Q s as in the original "bottom-up" paper [39]. Before we analyze the evolution of the background energy-momentum tensor in detail, we briefly review how the phase space distribution of gluons evolves in the "bottom-up" equilibration scenario [39], which anticipated the essential physics at weak coupling. Even at moderate values of the coupling constant (λ = 10), previous numerical simulations have shown that the system roughly follows the three different stages of the bottomup picture [37]. This is illustrated in Fig. 2, where we present contour plots of the phase space distribution as a function of time. At early times, shown in Fig. 2(a), the longitudinal expansion competes with the elastic 2 ↔ 2 broadening of the distribution function. Subsequently, soft bremsstrahlung emissions from high momentum partons begin to fill up the low momentum phase space as shown in Fig. 2(b). These soft particles can equilibrate via elastic 2 ↔ 2 interactions, leading to the formation of a soft and approximately thermal component at low momentum seen in Fig. 2(b); simultaneously the few remaining energetic partons ("mini-jets") are quenched and lose their energy to the soft thermal bath. Eventually, the mini-jets are fully quenched, the system approaches local thermal equilibrium, and the distribution function becomes well described by Bose-Einstein distribution with viscous corrections, c.f. Fig. 2(c).

C. Approach to hydrodynamics for the homogeneous background
Based on the above picture of the underlying microscopic dynamics, we will now discuss how the boost invariant homogeneous background evolves towards viscous hydrodynamics. As explained in Sec. II A, our main interest is in the non-equilibrium evolution of the energymomentum tensor, which can be calculated directly from the underlying particle distribution functionf (τ, p) where ν g = 2(N 2 c − 1) = 16 is the number of degrees of freedom for a pure glue plasma. We note that away from τ T id. /(4πη/s) 2nd hydro asympt. equilibrium, T µν does not satisfy hydrodynamic constitutive equations, and the energy-momentum tensor captures only the first moments of the non-equilibrium distribution function. However, the lowest moments of T µν provide sufficient information for the late time hydrodynamic evolution, and we will therefore neglect the effect of higher-order moments throughout this work.
Our results for the evolution of the diagonal elements of the background energy-momentum tensor 3 , are compactly summarized in Fig. 3, which shows the time evolution of the stress-energy tensor for different values of the 't Hooft coupling, λ = 10, 15, 20, 25. Motivated by previous work [37,46], we have rescaled the time and stress axes in order to fairly compare the physics at different values of the coupling. Specifically, at asymptotically late times the temperature approaches ideal hydrodynamics, parametrized as where Λ T is a dimensionful integration constant 3 Off-diagonal elements of the background energy-momentum tensor vanish by the symmetries of the underlying distribution. Hence off-diagonal contributions of T µν that may exist at initial time τEKT will be treated as linearized perturbations.
In Fig. 3 we first normalized the stress on the vertical axis by its asymptotic ideal hydrodynamics expectation and then rescaled the time on the horizontal axis by the equilibrium relaxation time τ R (τ ), which for typical modes is determined by the shear viscosity and the ideal temperature T id. (τ ; Λ T ) as After these rescalings the stress tensor follows a universal curve which is approximately independent of the coupling constant, at least for the range of couplings considered in this work. Such a scaling is guaranteed to work at late times where kinetic theory matches viscous hydrodynamics. Indeed, in second-order conformal hydrodynamics the energy density in a Bjorken expansion has the following asymptotic form (see Ref. [47] and Appendix A 1): where is a dimensionless combination of second order transport coefficients τ π , λ 1 and η/s. In leading-order kinetic theory all transport coefficients are functions of the coupling constant λ-the only free parameter of the theory. Consequently there is a one-to-one correspondence between the macroscopic parameter η/s and the microscopic parameter λ used in the simulation. (The details of the matching between the macroscopic and microscopic parameters are given in Appendix A.) Unlike η/s, the ratio C 2 has a much weaker dependence on the coupling constant [48], and thus the coupling constant dependence in the hydrodynamic result (Eq. (14)) is essentially completely contained in the scaled time variable 4 Returning to Fig. 3, we see that for a substantial range of coupling constants (or η/s), the non-equilibrium evolution of the energy-momentum tensor follows universal functions of scaled time τ T id. /(η/s), which are independent of the value of the coupling constant. Although such behavior is expected to emerge at late times, it is remarkable that the early time dynamics also exhibits such a universal behavior 5 . Similar observations have also been reported in [12,46,49], where this behavior is referred to as "hydrodynamic attractors". Exploiting the observed scaling property of the kinetic evolution, we will express the entire non-equilibrium energy-momentum tensor evolution in terms of a universal function of the scaled time τ T id. /(η/s). Since for a Bjorken expansion the longitudinal and transverse pressure can be readily determined from the conservation law P L (τ ) = −∂ τ (τ e(τ )) and the tracelessness of the energy-momentum tensor −e + 2P T + P L = 0, there is in fact only one independent function. Choosing the energy density e(τ ) as the independent variable, the entire pre-equilibrium evolution of the background energymomentum tensor can be parametrized as with the universal scaling function E x = τ T id. η/s . For asymptotically early times x = τ T id. η/s 1 the effective kinetic theory evolution is approximately fitted onto the free-streaming behavior, while for late times η/s is very well described by the hydrodynamic asymptotics, Eq. (14), with C 2 ≈ 1 determined in Appendix A. It is straightforward to construct a parametrization of E(x) which interpolates between the two limits and the resulting fit curve is displayed in Fig. 4. The explicit parametrization of E(x) is provided in Appendix A. We emphasize once again that, thanks to the scaling of the background evolution with η/s, the same fitted kinetic theory curve can be used for different values of η/s and different values of the initial energy density to map the early out-of-equilibrium energy density to the hydrodynamized energy-momentum tensor at later times.

D. Hydrodynamization time and pre-equilibrium entropy production
Based on the results in Fig. 3, we observe that, independently of the coupling constant, the kinetic description of equilibrating quark-gluon plasma overlaps well with hydrodynamics as soon as the reduced scaled time 5 We note that at very weak coupling this near equilibrium scaling ansatz will fail. Indeed, in the bottom-up scenario the hydrodynamization time scales parametrically as τ hydro ∼ α −13/5 s , whereas the hydrodynamic scaling ansatz predicts a parametrically larger time scale τ hydro ∼ α −3 s (c.f. Eqs. (17), (18) with η/s ∼ α −2 s ), due to the fact that the energy density of the hard modes is assumed to decay ∼ τ −4/3 rather than ∼ τ −1 . Of course, for moderate values of the coupling of practical interest the difference between α  is larger than unity τ T id. /(4πη/s) > 1. We therefore define a hydrodynamization time, τ hydro , as the boundary of applicability, i.e. τ hydro T id. /(4πη/s) = 1. Substituting the definition of T id. , we can express τ hydro in terms of η/s where Λ T is an energy scale defined through the asymptotic temperature Eq. (10), or alternatively to the asymptotic energy or entropy densities of the system. Specifically for a Bjorken expansion the temperature, energy density, and entropy density have the asymptotic forms: where Λ T , Λ E , and Λ S are all related to each other by the equation of state. We will parametrize all equations of state with an effective number of degrees of freedom where ν eff (T ) = ν g = 2(N 2 c − 1) = 16 for the gluon gas used in simulations, while for a three flavor gas of quarks and gluons ν eff (T ) = 47.5. Finally, at T ∼ 0.4 GeV lattice QCD simulations give ν eff ∼ 40 [50,51]. Using the definition Eq. (20) and thermodynamic identities, it is straightforward to show that the relation between the integration constants is 6 Hydrodynamic simulations usually adjust the energy density Λ E (or entropy density Λ S ) at equilibrium time to reproduce the multiplicity in the event.
which is consistent with other hydrodynamic simulations with a realistic equation of state (i.e. ν eff = 40) where τ s ≈ 4.1 GeV 2 [38]. Based on this estimate and using the hydrodynamization condition Eq. (18), we find that the hydrodynamization time of a boost invariant homogenous plasma is given by which provides a realistic bound for the applicability of relativistic viscous hydrodynamics (for a constant value of η/s = 2/(4π)). Changing the number of degrees of freedom from ν eff = 16 for a gluon gas to the more realistic ν eff = 40 increases the hydrodynamization time to 1.1 fm. 6 In Eq. (22) we used that (e + p)/e = 4/3, which is true within 5% even for lattice equation of state at T ∼ 0.4 GeV [50,51].
One additional consequence of the pre-equilibrium evolution is a rather rapid entropy production associated with the increase of the gluon number density per rapidity. With the full gluon distribution at our disposal we can immediately calculate the Boltzmann entropy s and the particle number n − (1 +f (τ, p)) ln(1 +f (τ, p)) , (25) Our results for the non-equilibrium entropy (s) and particle number (n) production are summarized in Fig. 5. We find that, independently of the coupling constant (or effective η/s), over 80% of total entropy per rapidity is produced by the end of the pre-equilibrium stage τ T id. /(4πη/s) ∼ 1. One observes that the nonequilibrium production can be very well reproduced for later times τ T id. /(4πη/s) 1 in second order hydrodynamics, provided a non-equilibrium entropy definition s non-eq (τ ) = s eq (τ ) (1 − τπ 4T η/s π µν π µν ) is used [47]. Finally, the entropy and gluon number densities roughly doubles from the start of kinetic evolution to the time τ hydro (see Fig. 5). Clearly, the rapid generation of entropy during the approach to equilibrium is very important in relating the initial state energy or gluon number density to the experimentally measured charged particle multiplicities (see Fig. 4 in Ref. [52] and Fig. 3 in Ref. [53], and references therein). Consequently, such large gluon multiplication factor needs to be taken into account in the correct estimation the properties of the initial state, e.g. the saturation scale Q s in Color Glass Condensate picture, from the measured multiplicities.

A. General decomposition of macroscopic response functions
Continuing the discussion of Sec. II A, we now look into the general properties of response functions for the linearized energy-momentum perturbations evolving on top of the out-of-equilibrium background. We consider only boost invariant perturbations in the transverse plane, and focus on the energy-momentum response to perturbations of the conserved charges-initial energy density δT τ τ and initial momentum density δT τ i .
By normalizing the perturbations to the background energy density T τ τ x (τ ), the evolution of energymomentum perturbations can be compactly summarized as Since we consider perturbations on top of a (locally) homogenous and boost invariant background, translation invariance guarantees that the response functions depend only on the difference x − x 0 and it is often more convenient to work in Fourier space, where we define the Fourier transformed response functionG µν αβ according to 7 Based on rotational symmetry in the transverse plane, one can further decompose the response functions into a tensor basis. For (scalar) energy perturbations the response functionG µν τ τ (k, τ, τ 0 , T τ τ x (τ 0 )) has only four independent structures while for (vector) momentum perturbations the decom-7 Note that vectors x and k are both confined to the transverse (x-y) plane.
position reads Since the longitudinal pressure componentsG ηη αβ are uniquely determined by the tracelessness of the energymomentum tensor, one is then left with a total of ten independent response functions, which need to be determined by a particular microscopic model.
Similarly to the discussion in k-space, the coordinate space response G µν αβ can be decomposed in tensors constructed from the radial vector r = x − x 0 . In practice, we first compute the response in k-space and then do the reverse Fourier transform [38] . The relations between momentum and coordinate space Green's functions are detailed in Appendix B.
B. Non-equilibrium response functions from effective kinetic theory The independent components of macroscopic response functions G µν αβ in Eqs. (29) and (30) need to be calculated by a particular microscopic theory. In this section we discuss the numerical realization of linear response in QCD kinetic theory around the non-equilibrium background presented in Sec. II C. At late times and close to thermal equilibrium, kinetic response functions are bound to approach the hydrodynamic limit, which is studied in Appendix D. Similarly, the early time dynamics can be profitably compared to the analytic results of collisionfree evolution, discussed in Appendix C.
We follow the methodology of Ref. [38] and linearize the phase-space distribution function f x,p around the backgroundf p , which is spatially homogeneous, but anisotropic in momentum space p = (p x , p y , p z ). We consider only boost invariant δf perturbations, which can be decomposed into a Fourier integral of plane wave perturbation δf k,p labelled by the wavenumber k in the transverse plane To linear order in perturbations, δf k,p evolves according to coupled Boltzmann equations where the collision kernel Since the evolution for different values of k on top of a homogenous backgroundf p decouples from each other [38], the linearized Boltzmann equation can be solved independently for each wavenumber k in the transverse plane. Even though the effective kinetic description of the pre-equilibrium dynamics requires the knowledge of the phase-space distribution δf k,p at the initial time τ 0 , one naturally expects the occurrence of memory loss during the evolution, so that the details of the initial phase-space distribution become irrelevant as the system approaches local thermal equilibrium. Since our ambition is merely to extract the energy-momentum tensor, a representative choice of δf k,p to characterize initial energy and momentum perturbations should be sufficient to describe the non-equilibrium evolution.
Based on a weak-coupling picture of the initial state, where the properties of the background distribution Hence one can motivate the initial phase-space distribution of scalar perturbations to be of the form where δQ s (k)/Q s = δT τ τ /T τ τ denotes the amplitude of the perturbation and the spectral shape is determined from the variation of the background distribution, Eq. (7), with respect to the scaleQ s . Similarly, considering that gradients of Q s (x) will lead to an initial velocity perturbation, e.g. v i (τ 0 ) ∝ −τ 0

∂iQs(x) Qs
, we can motivate vector perturbations of the form where v i (k) = δT τ i k,(v) /T τ τ denotes the amplitude of the initial velocity perturbation and the spectral shape is determined by a linearized velocity boost of the background momentum distribution. We note that in order to compute all response functions in the tensor decomposition in Eq. (30) it is important that we keep track of the independent components (labeled by the index i) of the momentum response. Even though at the level of the linearized Boltzmann equation, the actual magnitude of the perturbations is irrelevant in computing the response functions (as long as it remains sufficiently small to justify the linearized approximation at the relevant momentum scale), we find it convenient to choose an appropriate normalization. Defining the moments of the distribution function the corresponding energy and momentum perturbations associated with δf (s) k,p and respectively δf (v),i k,p at initial time are normalized such that for a highly oblate distribution (ξ 1) Given the above explicit form of the initial perturbations one can then determine the response of the energy-momentum tensor at any later time by numerically solving the kinetic equations, Eqs. (32) and (33). Once the solution to the linearized Boltzmann equation is calculated numerically, the response functions G µν τ τ (k, τ, τ 0 , T τ τ x (τ 0 )) can be directly constructed from the moments of the distribution function. For example, the energy response to energy or momentum perturbations is determined by the ratios Similarly, all the other components can be constructed by linear combinations of different components of δT µν as described in Appendix B

Scaling of response functions
We now present our numerical results for the nonequilibrium response functions calculated in effective kinetic theory. Even though generally the response func-tionsG µν αβ k, τ, τ 0 , T τ τ x (τ 0 ) depend separately on the wavenumber k, the initial and final times τ and τ 0 , the energy scale T τ τ x (τ 0 ), and the coupling constant The universal scaling functionG s s (see Eq. (42)) for the energy response to an initial energy perturbation as a function of the phase |k|(τ −τ0) at a fixed scaling time, τ T id. η/s . The different curves correspond to different coupling strengths (or η/s) which collapse onto a universal curve. The response functions at different scaling times τ T id. /(4πη/s) = 0.5, 1.5, 2.0 exhibit an equally good overlap (not shown). (b) The analogous plot for the energy response to an initial momentum perturbation, G v s . λ = g 2 N c , we expect that in analogy to the evolution of the background the number of independent variables can be drastically reduced by identifying appropriate scaling variables. Based on our analysis of the background evolution in Sec. II C, the natural candidate variables are the scaled evolution time τ T id. /(η/s) and phase |k|(τ − τ 0 ). Indeed we find that in the relevant range of parameters the postulated scaling property holds and the response functions can be compactly expressed in terms of a universal function of the scaling variables such that, for examplẽ Even though this behavior is expected to emerge in the hydrodynamic limit 8 of sufficiently late evolution times T id. τ /(η/s) 1 and for small wave-numbers |k|(τ − τ 0 ) 1, it is remarkable to observe that the scaling property holds across a much wider range of evolution times and wavelengths. As an illustrative example of the scaling, in Fig. 6 we present our results for the response functionsG s s andG s v given by Eqs. (40) and (41) at scaled evolution time τ T id. /(4πη/s) ≈ 1.0. Different curves in Fig. 6 correspond to simulations performed at different values of λ, which not only correspond to different effective η/s, but also amount to variations of the initial energy scale T τ τ x (τ 0 ), as seen from Eq. (7). Other components of the response function in the decomposition given by Eqs. (29) and (30) also show a good scaling with τ T id. /(η/s) and phase |k|(τ − τ 0 ) (not shown).
Because of the scaling of the response functions with η/s-the same kinetic theory response function computed for one set of initial conditions, can be used for different values of η/s and different values of the initial energy density to map the early out-of-equilibrium energy and momentum perturbations to the hydrodynamized energy-momentum tensor perturbations at later times. This is the procedure adopted in the pre-equilibrium propagator KøMPøST.
Finally, the coordinate space response functions used in the propagation formula Eq. (27) are obtained by Fourier transforming |k|-space components, e.g. Fig. 6, according to Eq. (28). The details of the procedure are given in Appendix B and also discussed in Ref. [38]. η/s and |x − x 0 |/(τ − τ 0 ).

Hydrodynamic constitutive equations
At late times when the system starts behaving hydrodynamically, the different components of response function decomposition Eqs. (29) and (30), are no longer independent, but are related by hydrodynamic constitutive equations. In second order hydrodynamics and for small wavenumber perturbations, the spatial part of energy momentum tensor δT ij can be written as a sum of energy δT τ τ and momentum δT τ i density perturbations, kinetic theory 1st order hydro 2nd order hydro kinetic theory 1st order hydro 2nd order hydro Comparison of the response functions for δT ij from an initial energy perturbation (seeG t,δ s andG t,k s in Eq. (29)) with the constitutive relations of first and second order hydrodynamics (Eqs. (43) and (44)).
with prefactors depending on first and second order hydrodynamic transport coefficients, i.e. η, τ π and λ 1 [38]. By comparing constitutive equations with the decomposition in Eq. (29), we find the following relation between response function components for initial energy perturbations In Fig. 7 we explicitly test the first and second order constitutive relations, Eqs. (43) and (44), in the kinetic evolution at scaled time τ T id. /(4πη/s) ≈ 2. Indeed for small wavenumbers |k|(τ − τ 0 ) < 6 the second order constitutive equations are well satisfied, demonstrating that the low wavenumber perturbations approach hydrodynamic regime at sufficiently late times τ T id. /(4πη/s) > 1. Similar constitutive relations can be derived for momentum response components and are given in Eq. (D19). The complete derivation is summarized in Appendix D.

IV. PRACTICAL IMPLEMENTATION: KØMPØST
Based on the general formalism and results of linearized effective kinetic theory presented in the previous sections, we will now describe a practical implementation of the pre-equilibrium evolution for hydrodynamic modeling of heavy-ion collisions-KøMPøST. Starting from a given profile of the energy-momentum tensor T µν (τ EKT , x) at an early time τ EKT , for example from the IP-Glasma model, we follow the procedure outlined below to calculate the energy momentum tensor T µν (τ hydro , x 0 ) at a later time τ hydro > τ EKT when the system is sufficiently close to local thermal equilibrium for viscous hydrodynamics to become applicable 9 .

a.
Decomposition in Background & Perturbations We first split the initial energy-momentum tensor T µν at KøMPøST initialisation time τ EKT in the background and perturbations where for linear evolution the decomposition into the background T µν x0 and perturbations δT µν is arbitrary, as long as the perturbations are sufficiently small. As discussed in Sec. II C we consider locally homogeneous boost invariant background energy-momentum tensor T µν x0 = diag (e, P T , P T , 1 τ 2 P L ) with only one independent component e(τ ). In order to obtain a smooth energy density profile from a discrete grid of input T τ τ values, we define the local background energy density as a Gaussian 9 Equilibration is not necessarily achieved everywhere at the same time τ , and the initial conditions could, in theory, be provided on a more complex (τ, x, y, η) hypersurface. However, it is the common practice to initialize hydrodynamic simulations on a constant τ hypersurface and we will follow this procedure in our present work.
weighted average around the point x 0 of interest where the Gaussian width is taken to be σ = ∆τ /2 10 . Here ∆τ = (τ hydro − τ EKT ) is the duration of kinetic evolution, i.e. the causal circle radius discussed in Sec. II A. Once the homogeneous diagonal background energy-momentum tensor T µν x0 is determined in the neighbourhood of point x 0 , the perturbation tensor δT µν x0 (τ, x) is obtained according to Eq. (45). In particular, the (small) initial off-diagonal components of energy-momentum tensor T µν are completely absorbed in the perturbation tensor, e.g. δT τ i x0 = T τ i . In accordance with the discussion in Sec. III A we only consider the response to initial energy δT τ τ x0 (τ EKT , x) and momentum perturbations δT τ i x0 (τ EKT , x). The initial perturbations in the shear-stress part of the energymomentum tensor, i.e. δT ij x0 are not taken into account in this work 11 .

b.
Background evolution & scale parameter Once the decomposition in background and perturbations is determined at each point x 0 in the transverse plane of the collision, we proceed to calculate the evolution of the background components of the energy-momentum tensor. Since in the relevant range of parameters the effective kinetic theory evolution exhibits a universal behavior in terms of the scaling variable τ T id.
η/s (see Sec. II C), we can immediately obtain the evolution for a specified values of η/s by matching the point associated with the initial energy density e(τ EKT ) = T τ τ x0 (τ EKT ) at time τ EKT to the universal scaling curve. Specifically, we determine the scale parameter Λ T (which fixes the temperature function T id. (τ ; Λ T ), see Eq. (10)) by solving the implicit equation where E(x) corresponds to the universal scaling curve for the evolution of the energy density (c.f. Sec. II C and Appendix A 2). Eq. (47) is the requirement that the initial time and energy density, i.e. τ EKT and e(τ EKT ), lie somewhere on the scaling curve shown in Fig. 4. Once the temperature function T id. (τ ; Λ T ) is known, the energy density at the hydrodynamic initialization time τ hydro can be read off from the same universal curve as Similarly, the background longitudinal and transverse pressure components, P L = τ 2 T ηη and P T = T ii , can be determined from the scaling curve as detailed in Appendix A 2, and we obtain the background energymomentum tensor T µν x0 at time τ hydro .

c.
Energy & momentum perturbations Next we propagate the initial energy and momentum perturbations, δT τ τ x0 (τ EKT , x) and δT τ i x0 (τ EKT , x), to calculate the contributions to all components of the full energymomentum tensor at hydrodynamic initialisation time τ hydro . We use the pre-calculated linear kinetic response functions G µν αβ discussed in Sec. III. First, the tabulated Fourier-space functionsG µν αβ τ T id.
η/s , |k|(τ − τ 0 ) are transformed to the coordinate space for the relevant values of scaled time τ T id.
η/s and radius |x − x 0 |/(τ hydro − τ EKT ), see Appendix B for details 12 . Subsequently, the contributions to the energy-momentum tensor at each point x 0 at the hydrodynamic initialisation surface τ hydro are determined by convoluting the coordinate space response functions with the initial energy and momentum perturbation as in Eq. (27). The coordinate space response functions contributing to δT µν (τ hydro , x 0 ) have only limited support, namely the neigbourhood of points x in the causal past of x 0 so in practice only a small number of spatial points contribute. We note that, according to the decompositions in Eqs. (29) and (30), we explicitly compute the contributions of energy and momentum perturbations to all components of energy-momentum tensor, without ever enforcing constitutive relations by hand. Adding the perturbations and the background produces the complete energy-momentum tensor at the end of KøMPøST evolution d. Decomposition of T µν in hydrodynamic variables Once the full energy-momentum tensor is obtained at hydrodynamic initialization time τ hydro , we perform a standard tensor decomposition into hydrodynamic variables [54,55]. Hydrodynamic simulations of heavy ion collisions describe the spacetime evolution of the energymomentum tensor T µν in terms of the energy density e, the flow velocity u µ , the bulk pressure Π and the shear-stress tensor π µν . The evolution of these fields is given by the energy-momentum conservation equation ∇ µ T µν = 0 and, in second-order Israel-Stewart formulations of relativistic viscous hydrodynamics [56], relaxation-type equations for the viscous fields, Π and π µν . The bulk pressure Π for conformal systems is exactly zero and the explicit expressions of T µν in terms of hydrodynamics fields is then given by where ∆ µν ≡ g µν +u µ u ν and the relation between energy and pressure is determined by the equation of state p = p(e) 13 . The tensor decomposition is done by identifying the local rest-frame velocity u µ and local energy density e as the time-like eigenvector and eigenvalue of T µ ν u ν = −eu µ . Once e and u µ values are known, π µν is obtained by tensor projection 14 . The independent hydrodynamic fields e, u µ and π µν are then passed to the subsequent hydrodynamic evolution.
Note that the linearized kinetic theory evolution does not guarantees the existence of a local fluid rest frame for arbitrary inputs; for sick cases, the procedure fails to find a meaningful rest frame. Such instances appear for the cases where the initial gradients are particularly steep (e.g. edges or peaks of the medium). Problems in extracting the flow velocity u µ from T µν in certain spatial regions are thus indicative of the linear approximation of the kinetic theory being pushed too far. Although these points make only a small fraction of the total points in the transverse extent of the fireball (quantified below), they tend to introduce instabilities in hydrodynamics code. Rather than attempting to address the problem in the hydrodynamic evolution, we developed a selective regulator of the kinetic theory output which we describe at the end of this section.
phase we use a lattice-based QCD equation of state [60], except when comparing to the conformal equation of state p = e/3. As for the first order transport coefficients, a constant shear viscosity over entropy density ratio η/s = 2/4π ≈ 0.16 is used, while the bulk viscosity is neglected. The second order transport coefficients are determined by relating them to the first order ones in the relaxation-time approximation [61,62]. The complete list of second order transport coefficients and hydrodynamic equations can be found in Ref. [63]. f. Hadronization At the end of hydrodynamic evolution the hadronic observables are computed from a constant temperature freeze-out surface with T FO = 145 MeV using the standard Cooper-Frye procedure [64]. We note that for simplicity only thermal hadronic observables, i.e. without hadronic decays, are used in this work, but the viscous corrections to the hadronic momentum distribution are taken into account as described in Refs. [59,63].
g. Regulators Below we document the regulator procedure for the limited instances when the energymomentum tensor T µν obtained at the end of the linear kinetic evolution cannot be inverted to define a local fluid restframe. The regulator we developed was motivated by free streaming, which is a meaningful and robust initial stage model. The regulator identifies regions of large gradients or low densities and drives the kinetic theory response towards free streaming in these regions by selectively lowering the scaling variable x = τ T id. (τ )/(η/s) to compute the kinetic response. We emphasize that while the regulator is important for making the hydrodynamic simulations run, it produces only minimal modifications of the hydrodynamic input and does not affect the physical results discussed in Sec. V. Thus, our pragmatic purpose here is to document the computer code.
Physically in heavy ion collisions the low density regions at the edges of the fireball can be never meaningfully described as a hydrodynamic medium. Hence there is no need to assume that η/s is constant throughout. For the purposes of estimating the scaling variable we therefore define (η/s)(T ) which grows large in regions of low density, i.e.
where C 0 = 100 MeV, and (η/s) 0 is a constant physical input parameter, typically of order 1/4π 15 . The scaling variable is then replaced by the number of relaxation times between times τ 1 and τ 2 Taking C 0 and τ 1 to zero, returns x 0 to the canonical scaling variable For obtaining the background energy density from the scaling curve, Eq. (17), we use the scaling variable value x 0 (τ hydro , 0), while for propagating the perturbations we use x 0 (τ hydro , τ EKT ), which provides a slightly better scaling parametrization of the Green functions 16 . Using x 0 as opposed to the canonical value, Eq. (54), removes a number of instabilities near the edge of the grid, but does not regulate the occasional regions of very high gradients in the central region of the fireball. The remaining instabilities arise when non-linearities become important. Examining when the T µν decomposition fails to find a rest frame, we determined (semiempirically) that for δT 00 the T µν decomposition may fail to find a rest frame. Certainly when z is of order 0.5 the linearized kinetic theory has reached its limit of applicability. For large z we regulated x 0 according to where S(z; z 1 , z 2 ) is a monotonic cubic spline interpolating between unity for z < z 1 , and zero for z > z 2 . In practice we take z 1 = 0.4 and z 2 = 0.7. When a regulated value x reg is used as opposed to x 0 the dynamics is pushed closer to free streaming limit in localised regions of steep gradients as shown in Fig. 8.
Technically the event is processed in two passes. In the first pass no regulator is used, and the scaling variable x 0 is calculated according to Eq. (53). From the T µν + δT µν of the first pass we record the size of the unregulated perturbations as measured by the variable z in Eq. (55), but we do not attempt to find a local fluid rest frame yet.
In the second pass we use z from the first pass to determine a regulated value of scaling variable x reg , Eq. (56). Then x reg is finally used to propagate the background and perturbations from τ EKT to τ hydro . The second pass T µν = T µν + δT µν has a rest frame decomposition which is passed on to the hydrodynamics code. We quantify the effect of the regulator by looking at the relative change in T τ τ component with and without a regulator 16 Typically τ hydro τEKT and the difference between the two values of the scaling variable is small. The relative change δ for MC-Glauber and IP-Glasma initial conditions for different evolution times is recorded in Table I. The transverse momentum flow grows with time and, according to criterion in Eq. (55), more points need to be regulated. However, from Table I it is clear that only a small fraction of points in the transverse extend of the fireball are affected and only for longest evolution times the change is at a few percent level.

V. EVENT-BY-EVENT PRE-EQUILIBRIUM DYNAMICS & MATCHING TO VISCOUS HYDRODYNAMICS
We will now illustrate the applicability of our framework to perform event-by-event simulations of the preequilibrium dynamics of high-energy heavy-ion collisions. Since the kinetic theory equilibration scenario described in the previous sections provides a smooth crossover from the early stage of heavy ion collisions to the viscous hydrodynamics regime, we will demonstrate with the example of two initial state models how initial conditions for hydrodynamic simulations can be obtained within our framework.
We first consider the Monte Carlo Glauber (MC-Glauber) model [17], which provides a phenomenological ansatz for the energy deposition in the transverse plane of heavy ion collisions, based on the location of binary nucleon collisions. Since MC-Glauber is a not a dynamical model, most phenomenological studies use an initialization time τ hydro ∼ 0.5 − 1 fm which is chosen empirically. However, we will show that the framework described in this work greatly reduces the sensitivity of the hydrodynamic evolution to the initialization time τ hydro .
In the second part of this section, we will also consider the IP-Glasma model [65,66], which provides a dynamical description of particle production and energy deposition. In this model color fields in each nucleus are sampled from a saturation model [67,68] and subsequently evolved with classical Yang-Mills evolution to times τ ∼ 1/Q s ∼ 0.1 fm. Since the early time dynamics of IP-Glasma matches smoothly onto our effective kinetic description, this implementation amounts to a complete dynamical evolution within a weak coupling framework.
While the microscopic IP-Glasma model provides an initialization for the entire energy-momentum tensor of the collision in 2+1D, the MC-Glauber model is typically used as an ansatz only for the transverse energy density 17 , without specifying the other components of the energy-momentum tensor. In the language of this paper, this means that IP-Glasma initial conditions provide both energy and momentum perturbations 18 , while the Glauber model only contains energy perturbations.
We note that the hydrodynamic initialization time τ hydro is treated as a variable in this section. The values of τ hydro used are of the order of the background hydrodynamization time given by Eq. 24, but not equal to it.

A. Energy perturbations with MC-Glauber initial conditions
We start our discussion with the Monte Carlo Glauber initial conditions, which provides an ansatz for the energy density distribution e(x) at each point in the transverse 17 It is also common to use Glauber model as an ansatz for the entropy density, which is then related to the energy density through the equation of state. In this work, the Glauber model is used as an ansatz for the energy density directly. 18 We note that IP-Glasma also provides higher order fluctuations, e.g. of the different T ij components. However, as discussed in Sec. III A we limit ourselves to the energy-momentum response to the fluctuations of conserved quantities like energy and momentum.
Even though the Glauber model itself does not provide an intrinsic time scale τ EKT for the dynamic evolution, it is clear that this time scale should be at least on the order of the formation time ∼ 1/Q required for semi-hard particles to go on-shell. Since we anticipate the typical momenta Q ∼ 2 GeV for central Pb+Pb events, we will use τ EKT = 0.1 fm in the following if not stated otherwise 20 . Similarly, the specific form of the energy-momentum tensor T µν = e(τ EKT , x) × diag (1, 1/2, 1/2, 0) in (τ, x, y, η) coordinates, can be motivated from the fact that due to the kinematics of high-energy collisions, the longitudinal 19 We used the publicly available T R ENTo code [18] to generate the event. In our study, the reduced nuclear thickness parameter p of the model was set to unity to obtain a typical participant Glauber model, while the negative binomial parameter for nucleon fluctuations was set to k = 1, and the nucleon smearing width w = 0.5 fm. 20 We checked explicitly that the sensitivity of our results to this choice is relatively small, as long as the initial energy density is rescaled by an appropriate factor, which can be deduced from the relation eτ = e 0 τ 0 + τ τ 0 dτ T zz for Bjorken expansion. Since at very early times T zz e, such rescaling effectively mimics a free-streaming evolution. momentum of each particle in the local rest frame is negligibly small compared to its transverse energy, such that the longitudinal pressure approximately vanishes at very early times (c.f. Sec. V B for a microscopic description of the early time dynamics). Since the T τ i momentum components of the energy-momentum tensor are also initialized as zero in the MC-Glauber model, the effective kinetic theory evolution of the energy-momentum tensor in Eq. (58) involves only energy perturbations.
We used a single MC-Glauber event with b = 0.9 fm impact parameter, N part = 408 participants and spatial eccentricity 2 = 0.064. The total transverse energy per unit rapidity in the event at τ EKT = 0.1 fm was set to d 2 x T τ τ = 5.3 × 10 4 GeV/fm to normalize the energy density distribution. The transverse distribution of a proxy quantity τ e 3/4 ∼ sτ for entropy per rapidity is shown in Fig. 9 for reference. Ultimately, after pre-equilibrium, hydrodynamic evolution, and freeze-out, this event corresponds to a midrapidity charged hadron multiplicity of dN ch /dη = 1870. The event is thus in line with a central LHC at √ s N N = 2.76 TeV, which, for reference, have dN ch /dη = 1601 ± 60 and N part = 383 ± 3 for the 5% most central events [69]. The energy momentum tensor in Eq. (58) features a large pressure anisotropy indicating that the system at τ EKT is still far from local equilibrium, and cannot be described properly by ordinary viscous hydrodynamics. However, the use of KøMPøST to describe the subsequent pre-equilibrium evolution (τ EKT < τ < τ hydro ) leads to the onset of hydrodynamic behavior that can be used as proper initial conditions for hydrodynamic evolution at τ hydro . The overlap in the range of validity of the pre-equilibrium and hydrodynamic phase ensures that the subsequent evolution is essentially independent on the switching time τ hydro . In practice, the smoothness of the transition from the early stage of heavy ion collisions to hydrodynamics can be quantified in multiple manners. In what follows, we look at averages and profiles of the hydrodynamics fields -or equivalently the energy-momentum tensor -as well as hadronic observables, and investigate their dependence on the hydrodynamic initialization time τ hydro .

Average hydrodynamic fields
Since realistic fluctuating initial conditions are used, hydrodynamic fields have a complicated profile in the transverse xy-plane. As a starting point, we look at transversely averaged values of hydrodynamic fields. We define averages . . . as where d 2 x denotes an integral over the transverse coordinates. The factor u τ is inserted for covariance, as would be obtained from the (covariant) surface flux dΣ µ u µ with dΣ µ in the (proper) time direction.
In Fig. 10(a) we present the evolution of τ e 3/4 , which is akin to the entropy per rapidity sτ of the boostinvariant system, as a function of physical time τ . Initial conditions at τ EKT = 0.1 fm are evolved with the kinetic theory in KøMPøST (c.f. Sec. IV), and subsequently passed to the hydrodynamic model at five values of τ hydro : 0.4, 0.6, 0.8, 1.0 and 1.2 fm. We first note that Fig. 10(a) shows the expected transition from τ e 3/4 being approximately constant at early time to subsequently dropping when the transverse expansion becomes important. More importantly one observes from Fig. 10(a) that the evolution of τ e 3/4 depends very weakly on the hydrodynamics initialization time τ hydro .
In Fig. 10(b) we show the rise of radial velocity The kinetic theory preequilibrium captures well the rapid rise in the radial flow at early times, which levels off to a steady radial increase at later times. It is again remarkable that the spread between the calculations for different values of τ hydro between 0.4 fm and 1.2 fm is at a percent level.
In order to follow the evolution of the azimuthal anisotropy, we computed the integrated transverse stress denotes an integral over the transverse plane without the energy weight, as T µν is already "energy weighted". Examining the principal axes of [T ij ] s , we define the momentum ellipticity which provides a measure of the elliptic flow as function of time.
Our results for momentum ellipticity are shown in Fig. 10(c), and indicate that the dependence of the averaged hydrodynamic fields on τ hydro is modest even if τ hydro is varied from 0.4 to 1.2 fm. Generally, the kinetic theory slightly under-predicts the pressure anisotropy ε p (τ ) in the hydrodynamic evolution. ε p (τ ) is inherently quadratic in the flow velocity, suggesting that including the first nonlinear couplings between the background radial flow and the generated elliptic flow could improve the agreement here.
The effective kinetic theory approach used in this work assumes that the system is conformal, i.e. the pressure is p(e) = 1 3 e. Even though QCD is nearly conformal at very high temperatures, deviations from conformality are expected for the range of temperatures of order 150−600 MeV encountered in heavy ion collisions at the RHIC and LHC, as can be observed from Fig. 11, where the ratio e/3p is shown for a QCD equation of state [60]. Since hydrodynamic simulations of heavy ion collisions necessarily require a realistic equation of state, this leads  to discontinuous matching of the energy-momentum tensor. Specifically, when the energy-momentum tensor of KøMPøST is decomposed and passed to the hydrodynamics, the energy-momentum tensor of the hydrody-namics becomes: which is not equal to T µν of the effective kinetic theory because p QCD (e) = p conformal (e). Unfortunately, there is no obvious way to improve on this procedure, as a better matching will ultimately require breaking the conformal symmetry in the pre-equilibrium phase which is of higher order in α s . On the other hand, it is straightforward to study and quantify the effects associated with this break of conformality, by replacing the QCD equation of state by a conformal one and reproducing Fig. 10(a-c). These results are presented in the bottom row of Fig. 10, where the different panels (d-f) again show the time evolution of τ e 3/4 , v ⊥ and ε p . It is clear from these figures that the τ hydro dependence, which was already small for a QCD equation of state, is even smaller with the conformal equation of state. In particular, all of the τ hydro dependence for τ e 3/4 -which amount to approximately 10% between τ hydro = 0.4 and 1.2 -is explained by the break of conformality between the initial conditions and the hydrodynamic evolution with a QCD equation of state. The flow observables v ⊥ and ε p also have a smaller dependence on τ hydro when a conformal equation of state is used, as can be observed by comparing Fig. 10 (b) with (e), and (c) with (f). The effect is not as significant as for τ e 3/4 , however, in part because the τ hydro dependence was already small in the first place with the QCD equation of state.

Transverse plane profiles of hydrodynamic fields
Based on our analysis in the previous section, the average energy, transverse velocity and momentum anisotropy were found to have a weak dependence on the value of τ hydro . The reason for this insensitivity to τ hydro is that after a short evolution, the dynamics of energy-momentum tensor in the effective kinetic theory approaches that of hydrodynamics, and the two descriptions are approximately equivalent. However, the integrated quantities shown in the previous subsection have a significantly reduced sensitivity to any type of fluctuations in the hydrodynamic fields. In this section we show profiles of hydrodynamic fields in the transverse plane of the collisions, to emphasize that the smaller features of the fields do not show any significant sensitivity to the hydrodynamic initialization time τ hydro either. Since the effect of the transition from a conformal kinetic theory to a non-conformal hydrodynamics evolution was quantified in the previous section, it is not revisited again here and all results that follow were obtained with a QCD equation of state.
In Fig. 12 we show profile plots of τ 3/4 (top row) and the flow velocity v x (bottom row) along the x-axis at midrapidity. Different panels (a-f) show the profiles at different hydrodynamic evolution times τ = 1.2, 2.0 and 5.0 fm. The different curves on each panel correspond to a hydrodynamic initialization times τ hydro from 0.4 to 1.2 fm. One observes that even for differential observables, the sensitivity to the initialization time of the hydrodynamic evolution is very small. Except for a few percent change in the overall normalization of τ 3/4 between the different curves, which can be attributed to the mismatch between the conformal equation of state in KøMPøST to the QCD equation of state in the hydrodynamic evolution, the profiles look essentially identical indicating a robust matching of the pre-equilibrium dynamics to viscous hydrodynamics.
One can further probe the approach of kinetic theory towards a hydrodynamic evolution by comparing the outof-equilibrium shear-stress tensor π µν from the kinetic theory evolution with an estimate from the Navier-Stokes value π µν = −ησ µν , where σ µν is calculated from the velocity profile, see Eq. (D7). In Fig. 13(a) we plot the value of π xx + π yy for τ hydro = 0.8, 1.0, 1.2 fm and Navier-Stokes estimate in dashed lines. We see that in most of the collision area the kinetic theory result approached hydrodynamic constitutive equations. One exception is the sharp edges of the fireball where the small gradient assumption breaks down. However, it is not clear that either hydrodynamics or linearized kinetic theoryà la KøMPøST provide an accurate description of the spacetime evolution of the edges. The regions where a good matching between kinetic theory and hydrodynamics is expected can be quantified with the typical momentum relaxation time τ R (τ ) defined at Eq. (13) in Section II C, using τ /τ R (τ ) 4π, i.e. τ T id. /(4πη/s) 1. This ratio can be calculated locally and indicates how the approach to hydrodynamics varies in the transverse plane. The result is shown in Fig. 13(b). As estimated in Sec. II D, we find that for starting times τ hydro > 0.8 fm, most of the medium is at scaled time τ T id. /(4πη/s) > 1 where hydrodynamics becomes applicable. Since the local energy density at the edges is significantly smaller, the edges of the fireball remain at τ T id. /(4πη/s) < 1 for a longer time, quantifying the statement that the approach to hydrodynamic behavior does not occur isochronously.

Hadronic observables
Based on the successful matching of the early time pre-equilibrium stage to the subsequent hydrodynamic regime discussed in the previous sections, we now investigate the impact of a consistent description of the early time dynamics on the final state hadronic observables computed after the freeze-out of the hydrodynamic evolution. We focus on the multiplicity dN π /dy, the average transverse momentum p π T and the v π 2 of thermal pions 21 , which can be thought of as analogues of the integrated hydrodynamic fields shown in Fig. 10. We note that, although only pion observables are shown in this section, we verified that similar results are found for heavier hadrons such as kaons and protons.
Starting with the same Glauber initial conditions at τ EKT = 0.1 fm, the effective kinetic theory is used to evolve the energy-momentum tensor up to five different times τ hydro from 0.4 fm to 1.2 fm, as in the previous sections. Subsequently, the hydrodynamic evolution is performed up to the isothermal freeze-out where hadronic observables are calculated. In Fig. 14, our results for the pion multiplicity dN π /dy, the mean transverse momentum p π T and the v π 2 are plotted as a function of τ hydro . In addition to the results obtained from the effective kinetic theory pre-equilibrium evolution, the dependence of hadronic observables on τ hydro is also shown for two other types of pre-equilibrium evolution: free-streaming 22 and simple Bjorken τ −4/3 scaling. The use of free-streaming to describe the pre-equilibrium dynamics has been studied previously in [21,22]. Similarly, the procedure of scaling the energy density of a set initial condition with τ −4/3 , as would be expected for a system undergoing ideal Bjorken hydrodynamic expansion, is also used regularly in heavy ion physics to rescale the energy density of the initial conditions when changing the initialization time of hydrodynamics.
In Fig. 14(a) we show pion multiplicity dN π /dy as a function of hydrodynamic initialization time τ hydro for the three different pre-equilibrium evolution scenarios. We find that for all pre-equilibrium scenarios, the multi- plicity has an approximately linear dependence on τ hydro . For the effective kinetic theory, the multiplicity is only approximately ∼ 5% smaller if the hydrodynamics is initialized at τ hydro = 1.2 fm rather than τ hydro = 0.4 fm, while in the case of the Bjorken τ 4/3 , this figure is significantly larger, ∼ 15%. Conversely, for a free-streaming pre-equilibrium dynamics the longitudinal pressure is underestimated during the pre-equilibrium phase, such that the energy decreases much less rapidly than in hydrodynamics, and the multiplicity is ∼ 8% larger with τ hydro = 1.2 fm than with τ hydro = 0.4 fm. We conclude that overall, the pion multiplicity has the smallest dependence on τ hydro when KøMPøST is used to describe the early stage of the medium evolution, although all curves are relatively flat.
In Fig. 14(b) we look at the radial flow dependence on τ hydro , again for kinetic theory, free streaming and τ 4/3 scaling. We find that for kinetic theory equilibration, the radial flow build-up is consistent with hydrodynamic evolution and the mean pion p π T is independent of switching time. In contrast, for the free streaming evolution, the overall energy scale grows slightly too rapidly, leading to a ∼ 3% change in p T if the hydrodynamics is initialized at τ hydro = 1.2 fm rather than τ hydro = 0.4 fm. Conversely, for Bjorken τ 4/3 scaling, no pre-flow is built up during the pre-equilibrium stage resulting in a decrease by ∼ 8% of p T for τ hydro = 1.2 fm rather than τ hydro = 0.4 fm. Similar observations can be made for the second order flow harmonic v 2 , presented in Fig. 14(c), albeit the overall magnitude of the variations is somewhat smaller in this case.
Besides the different τ hydro -dependence of p π T and v 2 , it is also clear from Fig. 14 that the p π T and v 2 ob-tained after a kinetic pre-equilibrium evolution is different from that obtained through the free-streaming and τ 4/3 scaling. Of course the same is also true for the multiplicity in Fig. 14(a): given the same initial conditions, dN π /dy is larger for a free-streaming evolution than a kinetic one, and smaller for τ 4/3 scaling. In practice, the normalization of the initial conditions of hydrodynamics simulations of heavy-ion collision is always adjusted so as to reproduce hadronic multiplicities. It is thus relevant to ask how p π T and v 2 change for free-streaming and τ 4/3 scaling if their respective initial condition normalizations are adjusted so that they produced the same pion multiplicity as the kinetic theory evolution. This result is shown with the open symbol on Fig. 14, for a fixed τ hydro = 1.2 fm. We see that even if the pion multiplicity is fixed by hand, different pre-equilibrium dynamics still lead to a difference in p T and v π 2 , although this difference becomes very small for free-streaming. In the case of Bjorken τ 4/3 scaling, the discrepancy remains relatively large, which we attribute in part of the lack of pre-equilibrium flow velocity resulting from this simplified scaling.

B. Energy and momentum perturbations with IP-Glasma initial conditions
Besides being applicable to general initial condition ansatzs, our framework of kinetic pre-equilibrium evolution can also be applied to microscopically motivated initial states. In this section we use the IP-Glasma model, where the evolution at early times of heavy ion collisions is described in terms of classical Yang-Mills (CYM) dy- namics [65,66]. We note that in contrast to our previous discussion of the MC-Glauber model, the microscopic IP-Glasma model also includes initial momentum fluctuations δT τ i , which are propagated by our kinetic theory evolution. Therefore we need the complete set of response functions discussed in Sec. III A.
Since classical-statistical field theory and effective kinetic theory have an overlapping range of validity, the combination of the two allows for a consistent weak coupling description of early time dynamics [37,44,70,71]. In principle, such combined approach describes particle production from the partonic structure of nuclei at high energies to the onset of hydrodynamics in the quarkgluon plasma.
In what follows, we use the original version of the IP-Glasma model, which is effectively 2+1D with boost-invariant fields in the longitudinal direction. Even though a qualitative matching between classicalstatistical field theory and effective kinetic theory has been demonstrated in previous works [31,37,39,72], no concrete implementation has been achieved to date for realistic full 2+1D simulation of heavy ion collisions. Because of the reduction to 2+1D (boost-invariant) fields, the dynamics in IP-Glasma becomes effectively freestreaming after τ ∼ 1/Q s ∼ O(0.1 fm) of classical Yang-Mills evolution [34]. This is different from the first stage of "bottom-up" equilibration, which is recovered in a full 3+1D classical-statistical simulation [31,45]. However in practice, we find that at weak coupling (larger values of η/s) the matching between classical Yang-Mills and kinetic theory is still rather smooth, as long as the switching to the kinetic description is performed at sufficiently early times (1/Q s τ EKT (4πη/s)/T id. ) as demonstrated below.
In order to illustrate the smooth matching at early times, we analyze the evolution of transversely averaged longitudinal P L and transverse P T pressures defined by Eq. (9), relative to the mean background energy density e in a particular event shown in Fig. 15. The IP-Glasma event used in our analysis is a central event with N p = 404 participant nucleons and an eccentricity 2 = 0.126 at time τ EKT = 0.2 fm. The total energy in the event is d 2 xτ EKT T τ τ = 4.4 × 10 4 GeV fm , corresponding to a charged hadron multiplicity of dN ch /dη = 1419. Fluctuations smaller than τ EKT have been smeared with a Gaussian.
In Fig. 16 we present the time dependence of the pres-sure to energy ratio P L,T / e during various stages of the evolution. In panel a) we see that at very early times (τ 1/Q s ) in the IP-Glasma evolution, the longitudinal pressure is negative due to the presence of strong longitudinal color fields [73]. However, on a time scale τ ∼ 1/Q s ∼ 0.1 fm the fields decay, and the longitudinal pressure approaches zero, except for a slight overshoot due to the residual pressure in the classical fields. At this time the 2+1D classical Yang-Mills dynamics becomes essentially free streaming, and one should then evolve the system with QCD kinetics in order to describe the subsequent approach towards equilibrium. Since at very early times, the expansion-dominated kinetic theory is also effectively free-streaming, the classical Yang-Mills and kinetic theory evolutions can be smoothly matched, provided τ EKT 1/Q s is small in units of the relaxation time, τ EKT T id. /(4πη/s) 1. Indeed, in Fig. 16(a) we see that IP-Glasma to kinetic-theory transition becomes increasingly smooth as 4πη/s at early times is increased from 2 to 8. The energy density is high at such early times, and this would naturally lead to larger values of η/s in a theory where the coupling runs. Of course the physics of the running coupling is absent in our leading order analysis where 4πη/s 2.
In Fig. 16(b) we see that during kinetic phase the pressure anisotropy decreases, and after the system is sufficiently close to local equilibrium τ T id. /(4πη/s) > 1, the subsequent evolution can be smoothly matched to second order viscous hydrodynamics, as discussed in Sec. II C. Fig. 16 represents the integrated stress and energies across a realistic IP-Glasma event, which approximately follows a 1D Bjorken expansion for the first τ 2.5 fm/c. After this time the 3D hydrodynamic expansion starts, and the (integrated) estimate P L / e 1 3 does not provide a useful measure of local thermal equilibrium in the flowing fluid.
Having investigated the different stages of the evolution, we will evaluate the transition between kinetic theory and hydrodynamics in greater detail by monitoring the stress tensor, and by checking that hadronic observables are independent of the crossover time τ hydro . The analysis parallels the MC-Glauber initial conditions described in Sec. V A, and the discussion will highlight the differences.

Average hydrodynamic fields
As in Sec. V A 1, we consider the transversely averaged hydrodynamic fields of energy τ e 3/4 , velocity v ⊥ and momentum eccentricity p (Eq. (61)) after the linearized kinetic pre-equilibrium evolution KøMPøST until τ hydro = 0.4, 0.6, 0.8, 1.0, 1.2 fm starting from IP-glasma initial conditions at τ EKT = 0.2 fm. The transversely averaged hydrodynamic fields are shown in Fig 17(a-c). The transverse average of τ e 3/4 and the transverse velocity v ⊥ are both showing an overall small dependence on τ hydro consistent with what was observed in Sec. V A 1 for MC-Glauber initial conditions. The momentum eccentricity p (Eq. (61)), on the other hand, is showing a larger dependence on τ hydro than in the Glauber case. We verified that similar dependence is obtained for p in peripheral IP-Glasma collisions with appreciable background ellipticity (not shown). The initial momentum perturbations for IP-Glasma initial conditions (which are absent in MC-Glauber initialization) are also propagated by the kinetic theory evolution, but they only make a minor contribution to averaged energy density τ e 3/4 and radial velocity v ⊥ at later times.
We note that both IP-Glasma and the kinetic theory are conformal, which means that the breaking of conformality discussed in Section V A 1 also occurs and increases the dependence on τ hydro for hydrodynamic evolution with realistic equation of state. In the bottom row of Fig. 17 we verified that using a conformal equation of state for the hydrodynamic phase slightly reduces the τ hydro dependence of both the transverse average of τ e 3/4 and of the transverse velocity v ⊥ , while leaving p relatively unchanged, as discussed in the previous section.
We also vary the transition time between IP-Glasma and the kinetic theory τ EKT , namely we do simulations with τ EKT = 0.1 fm and 0.2 fm, while keeping the crossover time to τ hydro = 0.8 fm fixed. As shown in Fig 18, all three averaged fields -the transverse energy, velocity and momentum eccentricities -show very little dependence on τ EKT . This is expected as at the crossover time the 2+1D Yang Mills evolution and kinetic theory are both close to free streaming. For the rest of this section we will use a fixed transition time of IP-Glasma to kinetic theory of τ EKT = 0.2 fm.

Transverse plane profiles of hydrodynamic fields
Next, we scrutinize the matching between KøMPøST with IP-Glasma initial conditions and hydrodynamic evolution by looking at the transverse energy and velocity profiles along the x-axis (y = 0) in Fig. 19. The upper row shows the transverse profile of energy density τ e 3/4 at different times in hydrodynamic evolution. Different lines represent a varying length of kinetic theory pre-equilibrium evolution with crossover times τ hydro = 0.4−1.2 fm. A good overlap of different curves indicates a smooth matching between kinetic theory and hydrodynamics. The small spread in energy can be attributed to the conformal breaking discussed in previous section, c.f. Fig. 17(a) and (d). In the bottom row of Fig. 19, we show the transverse velocity v x along x-axis. In the central region of the plasma we observe a smooth matching between kinetic theory and a full 2+1D relativistic hydrodynamics. For small energy densities at the edge of the medium (|x| 5 fm), the velocities are not as smoothly matched, but according to the discussion in Sec. II D, these regions do not satisfy the criterion of hydrodynamization anyway. Although at later times the pre-equilibrium flow is dominated by the response to ini- tial energy gradients, the IP-Glasma initial conditions have a non-zero initial velocity at τ EKT = 0.2 fm, which contribute to the fine details of profiles shown in Fig. 19.
To test the convergence of viscous components of energy momentum tensor to hydrodynamic expectations, in Fig. 20(a) we show the shear stress tensor π µν profile for three values of the hydrodynamic initialization time τ hydro = 0.8 fm, 1.0 fm and 1.2 fm, and compare to the estimated Navier-Stokes value ∼ ησ µν (obtained from the velocity profile). The agreement with the Navier-Stokes value is not as good as observed for the smoother Glauber initial conditions (Fig. 13), although still reasonable. In panel (b) we show the local scaled evolution time τ T id.
η/s for different crossover times, which shows that the central part of the collision is within the hydrodynamic regime τ T id. /(4πη/s) > 1 at time τ hydro > 0.8 fm as stipulated in Sec. II D.

Hadronic observables
After checking the smooth matching between individual stages of the evolution, we now test the effect of the hydrodynamic initialization time τ hydro on the final state observables. To recap, the IP-Glasma initial conditions are evolved until τ EKT = 0.2 fm at which point the background energy density together with energy and momentum perturbations are passed to KøMPøST. After linearized kinetic theory pre-evolution the hydrodynamic fields like energy e, velocity u µ and shear-stress tensor π µν are passed to viscous hydrodynamic simulation at time τ hydro = 0.4−1.2 fm. Then thermal hadronic observables are computed at constant temperature T FO = 145 MeV freeze-out surface via the standard Cooper-Frye procedure. In Fig. 21 we show the thermal pion multiplicity dN π /y, the mean transverse momentum p π T and v π 2 , as a function of the hydrodynamic initialization time τ hydro . For comparison, we replace the kinetic theory evolution with free-streaming background evolution and response functions. As observed in the Glauber case (see Sec. V A 3 and Fig. 14), hadronic observables show very little dependence on τ hydro when a kinetic theory pre-equilibrium evolution is used. The pion multiplicity changes by less than 4%; for free-streaming preequilibrium the change is twice as large in the opposite direction. The mean radial p T is essentially independent of τ hydro , but is slightly increasing for free streaming pre-equilibrium. Finally, the elliptic flow slightly changes for both types of pre-equilibrium evolutions, but it is still smaller for kinetic theory description. All in all, the dependence on τ hydro is small for the KøMPøST preequilibrium evolution and the dependence increase when the kinetic theory is replaced by free-streaming description. We highlight that the τ hydro dependence of the pion v 2 is still small, despite the somewhat larger dependence of p seen in Fig. 17.

VI. EFFECTIVE DESCRIPTIONS OF EARLY TIME DYNAMICS
So far we have demonstrated the practical performance of our framework to describe early time dynamics of highenergy collisions. We will now investigate in more detail theoretical relations and practical comparison to other approaches previously discussed in the literature [15,[21][22][23][24].

A. Long wavelength limit of kinetic theory response
In viscous hydrodynamics small scale fluctuations dampen rapidly, and many final state observables are only sensitive to long wavelength perturbations [74][75][76][77][78][79][80]. Consequently, one could expect that also the preequilibration discussed in previous sections could be well captured by the long wavelength response. We will now discuss how to formalize and test this idea, based on a low |k| expansion of our non-equilibrium linear response formalism. Interestingly, it will also be useful to establish the relation of our formalism to previous ideas related to the concept of a "universal pre-flow" [23]. Details of the derivations presented in this section are worked out in Appendix E.
In order to study the low |k| limit of kinetic theory preequilibrium evolution to initial conditions of heavy ion collisions, we first filter out the small wavelength perturbations by performing a Gaussian smearing of the initial energy-momentum tensor T µν (τ 0 , x 0 ). Specifically, we define a coarse grained energy-momentum tensor with the same Gaussian smearing kernel S σ (x 0 − x 0 ) used to define the local background energyT µν x (τ 0 ) (c.f. Eq. (46)). Considering a space-time point (τ, x) on the future hydrodynamic surface, we can then perform our usual decomposition of the energy-momentum tensor into local backgroundT µν x (τ ) and perturbations δT µν x (τ 0 , x 0 ). However, instead of considering fluctuations on all scales, the initial energy-momentum tensor perturbation δT µν x (τ 0 , x 0 ) can now be decomposed further into short and long-wavelength components as where, recalling that the definition of the local back-groundT µν x (τ 0 ) involves the same coarse graining procedure, the long wave-length components of δT µν x (τ 0 , x 0 ) are then entirely given in terms of smooth fieldsT µν . If the initial profile of the energy-momentum tensor T µν (τ, x) is sufficiently smooth on length scales smaller than σ ∼ |τ − τ 0 | the short wave-length component is very small and can safely be ignored 23 . Neglecting short wavelength components in the following, the energymomentum tensor at later times is then given by 23 Note that even if the initial energy-momentum tensor features fluctuations smaller on length scales smaller than σ, short wave length components on scales less than σ visc ∼ (τ − τ 0 )/ τ T id. η/s , are subject to strong viscous damping and could still be neglected.
long wavelength coefficients 22. Linear and quadratic long wavelength response coefficients to initial energy perturbations (a) and initial momentum perturbations (b) from non-equilibrium kinetic theory evolution. Dashed lines show comparison to viscous hydrodynamic asymptotic derived in Appendix D 2.
By expressing the convolution in Fourier space, and expanding the response functions in powers of the wave number |k|, as and similarly for the other components, it is then straightforward to show that the long wave length response in Eq. (65) is proportional to the gradients of the smoothened energy-momentum tensorT µν (τ, x). Saving the details of the derivation for Appendix E, the energymomentum response to initial energy gradients is then given by which is a generalization of the previous result for the transverse pre-flow derived in [23,38]. Eq. (67) says that in addition to the "pre-flow" which develops during the equilibration process due to gradients [23], local maxima of the initial energy density, which have negative second derivatives, are depleted during the evolution. We can also obtain the long wavelength response to initial mo-mentum perturbations We extract the 1st and 2nd order coefficients in |k|(τ − τ 0 ) by performing a polynomial fit to the first few |k| values of kinetic theory response functions, e.g. Fig. 6. Our results for the zero momentum response G 0 s/v as well as the long wave-length coefficients s (n) ,v (n) from effective kinetic theory simulations, are compactly summarized in Fig. 22, where we plot the various response coefficients as a function of the scaling variables τ T id.
η/s . We find that the time dependence of the coefficients is rather weak, such that in practice the long wave length response can be approximated rather well by the hydrodynamic asymptotics, see Appendix D 2.
Before we compare the long wavelength results to the full treatment in KøMPøST, it is important to point out that the above separation into long and short wavelength modes introduces an artificial regulator dependence not present the full treatment. Specifically for the long wavelength filterS σ (k) ∝ e − σ 2 2 k 2 1 − σ 2 2 k 2 in Eq. (63), the first difference formally appears at quadratic order O(k 2 ). However, the O(k 2 ) can always be absorbed into an additive renormalization of the long wavelength coefficients, i.e. by subtracting σ 2 /(τ −τ 0 ) 2 from the quadratic s (2) s coefficient in energy response, as explained in Appendix E 2. Based on this renormalization scheme, the regulator dependence enters only at cubic order O(|k| 3 ) and the residual dependence can be used to quantify the systematic uncertainties of the long wavelength approximation.

B. Comparison of effective kinetic theory, long wavelength response & free streaming
We now turn to the comparison between the full kinetic theory and the low-|k| limit described in previous section. Our results are summarized in the top row of Fig. 23, where we compare profiles of the energymomentum tensor at the end of the pre-equilibrium evolution at τ hydro = 1.2 fm initialized with the same MC-Glauber initial conditions at τ EKT = 0.1 fm. We assess the robustness of the long wavelength results by using two different smearing widths σ 1 = σ 0 and σ 2 = 2σ 0 , where σ 0 = (τ hydro − τ EKT )/2 is the same Gaussian width used to define the local averaged background in the full kinetic theory response.
As it is visible from Fig. 23(a) the simplified low-k evolution accurate to quadratic order in small |k| reproduces the energy density profile rather well and is largely insensitive to the smearing width. However, most of the energy is evolved as a background according to universal scaling curves, Eq. (17), which is the the same for both the full kinetic and low |k| response.
In Fig. 23(b) we compare the transverse velocity profile in kinetic theory and long wavelength limit. Formally, the low-k limit given by Eq. (67b) is only accurate to linear order in k (c.f. [23]). However for the particular choice of the regulator σ 1 = σ 0 , the long wavelength limit approximately reproduces the cubic order flow response (see Appendix E 2 for details) and the resulting curve is very close to the full kinetic theory result. Despite this accidental agreement, it is also evident from the comparison with the curve for σ 2 = 2σ 0 that the long wavelength result exhibits a strong regulator dependence, which points to the fact that the actual leading order (∼ k) velocity profile provides a less accurate approximation of the full kinetic theory result.
Similar features can be observed in Fig. 23(c), where we compare the kinetic theory evolution of the shearstress tensor π µν with the corresponding long wavelength result obtained for simplicity through hydrodynamic constitutive equations. At leading viscous order π µν ∝ ησ µν is mainly dominated by velocity gradients, therefore the agreement is better for the low-|k| response with accidental cubic accuracy in velocities, but is also reasonably good for the different choice of coarse-graining.
It is constructive to compare in the same fashion the kinetic theory response with a simple model of the preequilibrium evolution prescription based on free streaming [21,22]. For completeness the free-streaming response functions are summarized in Appendix C. The comparison is shown in Fig. 23(d-f). In free streaming evolution the longitudinal pressure is completely neglected and thus the energy density decreases more slowly as a function of time, overshooting the kinetic theory results (see Fig. 23(d)). For the special case of initial (scalar) energy perturbations, the free streaming response functions for transverse velocity have approximately the same low-|k| expansion as the full kinetic theory evolution, see Eq. (D57). Therefore the transverse velocity profile in Fig. 23(e) between kinetic theory and free streaming response is almost indistinguishable. However, the physical momentum is not correct, because of the incorrect background energy evolution, (c.f. Fig. 23(d)).
Finally in Fig. 23(f) we compare the shear-stress tensor π µν in free streaming (the dashed-dotted line) to the full kinetic theory result (the solid line). In free streaming the longitudinal pressure is completely neglected, and thus the transverse stress π xx +π yy is too large when the system is passed to hydrodynamics. (We have reduced the free π xx +π yy by a factor of two in Fig. 23 for visibility.) If instead of the free streaming result for π xx +π yy , one just uses the free energy and velocity profiles, Fig. 23(d-e), together with the hydrodynamic constitutive relations, the velocity v x (center) and shear stress (π xx + π yy ) (right).
shear-stress tensor is closer to the kinetic theory result, as indicated by a Navier-Stokes estimate ∼ ησ µν in the figure.
To summarize, the low-|k| limit does a good job at describing the velocity, but it makes no predictions for the background energy density as a function of time. Therefore it must be supplemented by a theory (such as QCD kinetics) which describes the longitudinal pressure at early times. Free streaming also correctly describes the velocity response, but in this case the energy density and second order hydrodynamic variables must be continually readjusted in order to have a smooth transition to hydrodynamics. On a practical note, using KøMPøST is computationally no more expensive than the other alternatives, and it should become the pre-hydro engine of choice.

VII. SUMMARY & OUTLOOK
In this paper we developed a linear response framework to describe the non-equilibrium evolution of the energy-momentum tensor during the early stages of highenergy heavy-ion collisions. Based on microscopic input from effective kinetic theory simulations, we presented a practical implementation of the "bottom-up" thermalization scenario with realistic fluctuating heavy-ion initial conditions and demonstrated a consistent matching between the pre-equilibrium and the hydrodynamic evolution on an event-by-event basis. Our linear kinetic preequilibrium propagator KøMPøST [41] provides a practical implementation of a systematic procedure to propagate initial out-of-equilibrium perturbations to the hydrodynamic initialization time [38]. Crucially, for short evolution times the equilibration dynamics can be described as sum of local background equilibration and linear response to fluctuations of initial conserved energy and momentum density, Eq. (4). We developed a concrete realization of the general linear response formalism presented in Sec. II A and Sec. III A using QCD kinetic theory with gluonic degrees of freedom and extrapolation to moderate values of the coupling constant λ [30,37,38]. Within our framework, the entire kinetic equilibration of T µν is compactly summarized by a single evolution curve of the homogeneous background (c.f. Fig. 4) and a handful of response functions (c.f. Figs. 26 and 27 ). Interestingly, we find that for the relevant range of coupling constants, corresponding to realistic values of η/s, the kinetic equilibration only depends on the scaling time ∼ τ T /(η/s) such that pre-tabulated kinetic response functions can be used for event-by-event simulations.
We found that for typical QGP parameters the leading order kinetic theory predicts a hydrodynamization time around τ hydro ∼ 1 fm (c.f. Eq. (24) and Ref. [38]). During this time the initial gluon number density per rapidity roughly doubles and, thus, significantly alters the relation between final particle multiplicities and the entropy density per rapidity in the initial state (see Fig. 5).
We applied our formalism to two widely used initial state descriptions: a phenomenological MC-Glauber ansatz for transverse energy density deposition in Sec. V A and the first-principle dynamical IP-Glasma model (which contains both energy and momentum fluctuations) in Sec. V B. In both cases we demonstrated a smooth matching between the kinetic theory and hydrodynamics, both in the hydrodynamics fields, e.g. Fig. 12, and in the final hadronic observables, see Figs. 14 and 21, largely eliminating the dependence on the crossover time τ hydro . Note that the kinetic response functions reproduce the spatial components of the energy-momentum tensor δT ij without explicitly imposing the constitutive relations, (c.f. Figs. 13(a) and 20(a)).
Finally we studied the kinetic response in the "low wavelength" (hydrodynamic) and "fast expansion" (freestreaming) limits to establish connections with the existing literature [15,[21][22][23][24]. We find that the first few terms in small |k|-expansion are sufficient to capture most of kinetic response (Fig. 23). Incidentally, the leading order velocity response in free-streaming limit agrees with kinetic theory, but energy density and shear-stress tensor evolution are described incorrectly.
Our publicly available kinetic propagator package KøMPøST [41] offers a simple, yet non-trivial description of hydrodynamization in heavy ion collisions. The reduced sensitivity on the hydrodynamic initialization time τ hydro brings us closer to ab initio initial conditions of heavy ion collisions. Moreover, there are multiple possibilities to build on this work to achieve an even better description of the early stage of heavy ion collisions. The inclusion of quarks as degrees of freedom would lead to a more realistic description of the pre-equilibration stage of collisions, which should provide a better matching with the QCD equation of state. Studying non-trivial back-ground evolution in kinetic theory, e.g. radial gradients, could improve the applicability of presented framework at the edges of the QGP fireball and in small systems.
More importantly, the formalism derived in this work to propagate linear perturbations on top of a smoother background can be used with response functions computed in limits other than weakly-coupled effective QCD kinetic theory. By using this framework to compare systematically the macroscopic description of equilibration from a weakly-coupled regime and a strongly-coupled one, one can hope to better constrain the real dynamics of the medium produced in heavy ion collisions.

ACKNOWLEDGMENTS
The authors would like to thank Björn Schenke for insightful discussions and for his help adapting the hydrodynamics code MUSIC for this work, and Liam Keegan for his contributions at the beginning of this project. Useful discussions with Jürgen Berges, Stefan Flörchinger, Yacine Mehtar-Tani, Klaus Reygers, Raju Venugopalan are gratefully acknowledged. Results in this paper were obtained using the high-performance computing system at the Institute for Advanced Computational Science at Stony Brook University. This work In this section we summarise the procedure of extracting the transport coefficients in kinetic theory for different values of the coupling constant λ. At asymptotically weak coupling, the transport coefficients can be obtained from perturbative calculations [48,81], however in this work we consider moderate values of λ. So instead of using perturbative formulas, we extract the transport coefficients directly from the evolution of homogeneous and longitudinally expanding plasma. Effectively we treat λ as a dial parameter in our leading order kinetic theory description to tune to different values of transport coefficients (which, perhaps, correspond to different value of the coupling constant in a higher order treatment).
First, the transport coefficients are obtained from the relaxation of pressure anisotropy. For Bjorken expansion and in second order conformal hydrodynamics one can show that it takes the following form (see Ref. [47] and Appendix D) where T (e) is the equilibrium temperature obtained from the Landau matching condition, Eq. (20), and C 2 is a constant combination of first and second order transport coefficients, see Eq. (15). The shear over entropy ratio η/s and C 2 are then fitted to the late time tail of pressure anisotropy. The fit results for different values of λ is shown in Fig. 24. Note that C 2 is only modestly dependent on the coupling constant and rather close to the weak coupling estimate C 2 ≈ 1.0 (τ π ∼ 5.1η/sT , and λ 1 ∼ 0.8ητ π [48]), so that in practice the same calculated value C 2 ≈ 1.0 can be used for the entire range of the coupling constants.
As discussed in Sec. II C, the hydrodynamic evolution collapses to the same scaling curve if time is counted in relaxation time units, see Eq. (13). In order to define a simple single parameter temperature scale at all times we introduced temperature T id. (τ ; Λ T ) in Eq. (10). Relating T id. to the equilibrium temperature T (e) at first viscous order Λ 2/3 T ≡ T id. τ 1/3 = τ 1/3 T (e) 1 + 2 3 η/s τ T (e) (A2) gives a simple and efficient way of determining the energy scale Λ T . Note that the pressure equilibration, Eq. (A1), can be now written explicitly as a function of time and three fit parameters: η/s, C 2 and Λ T . Remarkably, the pressure anisotropy evolution in Fig. 24 collapses to the same curve even at earlier times than when the second order asymptotics becomes applicable.
Finally for completeness we list the asymptotic expansions for temperature, energy and entropy at second order in terms of scaling variable τ T id. η/s : η/s τ T id. − 2 9 C 2 η/s τ T id.

Parametrization of background evolution
Below we provide explicit parametrizations of the evolution of the background components of the energy momentum tensor. Based on our discussion in Sec. II C, the kinetic theory evolution can be parametrized in terms of a single scaling function E(x) for the energy density We parametrize the scaling function E(x) by the following ansatz where step(x) is a smooth function interpolating between step(x 1) = 1 and step(x 1) = 0. Explicitly we take step(x) ≡ 1 2 1 + tanh At late times we recover the hydrodynamic limit Eq. (A5) which is quite well described by C 2 = 1.02 [48]. At early times, the energy evolution is matched to a freestreaming behavior The remaining fit parameters are √ 4πF 0 = 1.0767, 4πb = −0.284, (4π) 2 c = 0.134 (A12) while the step function is parametrized by x Switch /(4π) = 0.65, x Range /(4π) = 0.25 (A13) The transverse and longitudinal pressure is determined by functions P L (τ ) = ν g π 2 90 T 4 id. P(x), P T (τ ) = ν g π 2 90 T 4 id.
Here the auxiliary function P(x) is determined by the equation of motion for Bjorken expansion, i.e. ∂ τ (τ e) = −P L , A comparison of the parametrization with the simulation data is compiled in Fig. 25. We note that since the scaling of kinetic evolution curves is not perfect, especially for the longitudinal pressure at very early times, there is a systematic error in this parametrization at the order of a few percent. Similarly, for (vector) momentum perturbations, the decomposition take the form δ v (τ, τ 0 , |r|)+ 1 2 r i |r| δ jk + r j |r| δ ik G m,r v (τ, τ 0 , |r|) + r i r j r k |r| 3 G t,r v (τ, τ 0 , |r|) , with the different components obtained from Eq. (30) according to following transformations G s v (τ, τ 0 , |r|) = 1 2π d|k||k|J 1 (|k||r|)G s v (τ, τ 0 , k) , tum perturbations. We employ co-moving coordinates denotes as with i = 1, 2 labeling the transverse Lorentz indices. Our convention for the metric is mostly positive, i.e.

Second order viscous constitutive relations
Hydrodynamic constitutive relations take the form T µν = (e + p) u µ u ν + p g µν + π µν , where at second order the shear stress tensor π µν is given by 25 π µν = −ησ µν + ητ π u λ ∇ λ σ µν + 1 3 σ µν ∇ λ u λ where ∇ α u β = ∂ α u β − Γ µ αβ u µ denotes the covariant derivative and ∆ µν = g µν + u µ u ν is the usual projector to the rest frame. Based on these definitions, the components of the shear-stress tensor can be evaluated according to the following relations (no summation over