Matching the Nonequilibrium Initial Stage of Heavy Ion Collisions to Hydrodynamics with QCD Kinetic Theory

High-energy nuclear collisions produce a nonequilibrium plasma of quarks and gluons which thermalizes and exhibits hydrodynamic flow. There are currently no practical frameworks to connect the early particle production in classical field simulations to the subsequent hydrodynamic evolution. We build such a framework using nonequilibrium Green's functions, calculated in QCD kinetic theory, to propagate the initial energy-momentum tensor to the hydrodynamic phase. We demonstrate that this approach can be easily incorporated into existing hydrodynamic simulations, leading to stronger constraints on the energy density at early times and the transport properties of the QCD medium. Based on (conformal) scaling properties of the Green's functions, we further obtain pragmatic bounds for the applicability of hydrodynamics in nuclear collisions.

Collisions of heavy nuclei at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) heat up nuclear matter sufficiently to produce a plasma of deconfined colored degrees of freedom -the quarkgluon plasma (QGP) [1][2][3][4]. The properties of this deconfined plasma can only be constrained indirectly from the mass and momentum distribution of the final shower of color-neutral particles reaching the detectors. Extensive comparisons with measurements indicate that the spacetime evolution of the QGP can be described with relativistic viscous hydrodynamics [5][6][7][8][9] starting from a time τ ≡ τ hydro ∼ 1 fm/c after the collision, with other models describing the dynamics of the collision before and after this hydrodynamic phase. The success of hydrodynamic models has made it possible to study certain finite temperature properties of quantum chromodynamics (QCD), such as the specific shear viscosity η/s. Most remarkably it was found that η/s is on the order of a tenth in units of /k B (with k B being the Boltzmann constant), perhaps the smallest specific shear viscosity ever measured. Because constraints on the plasma's properties are obtained from multistage simulations, an intricate interdependency exists between the hydrodynamic evolution and the description of the earlier phase of heavy ion collisions. Obtaining more precise constraints on η/s and other transport coefficients of QCD hinges to a considerable extent on attaining a better understanding of this early stage of the collisions and its transition to hydrodynamics.
Before the time τ hydro , the deconfined matter is not amenable to a coarse-grained description in terms of macroscopic hydrodynamics fields, and a microscopic description must be used. Modeling the earliest phase of heavy ion collisions and its transition to hydrodynamics from first-principles QCD remains a formidable challenge. In this Letter, we address this challenge by showing how an effective kinetic theory (EKT) of weakly coupled QCD [10] can be used to smoothly describe the evolution of a general out-of-equilibrium energy-momentum tensor specified at very early time τ ekt τ hydro to its late time hydrodynamic form.
Naturally, our weakly coupled approach is particularly well suited for collisions of large nuclei at high energies, where a quantitative theory of gluon production can be developed based on the color glass condensate theory of high-energy QCD [11][12][13][14]. Based on a separation into slow and fast degrees of freedom, the early time dynamics of the collision in the color glass condensate approach is described in terms of classical gluon fields. Classical Yang-Mills equations determine the dynamics of the sys-tem until a time τ ekt at which the phase space density of gluons (per ) becomes less than ∼ 1/α s (Q s ), with α s being the strong coupling constant evaluated at the saturation scale Q s . Subsequently, for τ > τ ekt kinetic processes dominate the spacetime evolution of the highly occupied and highly anisotropic plasma of gluons [15,16]; detailed simulations of the classical Yang-Mills dynamics [15][16][17] confirm the onset of the "bottom-up" thermalization scenario developed in a seminal paper many years ago [18]. Since then there have been a multitude of analytical and numerical works devoted to clarifying the non-Abelian field dynamics in the earliest stages of the evolution [15][16][17][19][20][21][22][23][24], and to developing tools to simulate the subsequent approach towards local thermal equilibrium using QCD kinetics [25][26][27][28][29][30][31], involving several unique features such as non-Abelian collinear radiation [32] and dynamical screening [10]. In spite of this theoretical progress, a practical tool to bridge between the early time dynamics of strong color fields and successful hydrodynamic simulations at late times has not been achieved so far.
Based on a nonequilibrium linear response formalism [31] developed in detail in our companion paper [33], our Letter provides for the first time a concrete realization of a satisfactory theoretical description of the early time out-of-equilibrium dynamics. The computer code, called KøMPøST after its authors, is publicly available [34]. While the underlying kinetic approach can be justified rigorously for the collision of large nuclei only in the limit of very weak coupling, we show that the kinetic response to a variety of initial conditions can be smoothly extrapolated to physically relevant couplings by using an appropriate scaling variable. This makes KøMPøST a practical tool for RHIC and LHC phenomenology. Nonequilibrium linear response theory.-Using a nonequilibrium linear response formalism, the energymomentum tensor of the system is evolved from its nonequilibrium form at τ ekt up to a time τ hydro when hydrodynamics becomes applicable -see Fig. 1. Our approximation scheme is based on the separation of scales seen in Fig. 1 [31]: to determine the energy density at a spacetime point (τ hydro ,x), one needs only to propagate the nonequilibrium initial condition from points (τ ekt ,x ) in the causal past |x − x | < c(τ hydro − τ ekt ). Since this causal circle is small compared to the system size ∼ 2R, the energy-momentum tensor in this domain can be divided into a local average and (small) perturbations [35] In practice, the background T µν x (τ ) is calculated by a spatial average of T µν (τ ekt , x ) over the causal circle, and is assumed to be boost invariant and locally homogeneous in the transverse (xy-)plane. Differences between the full energy-momentum tensor T µν (τ ekt , x ) and the background T µν x (τ ekt ) are treated as linearized perturba- tions δT µν x (τ ekt , x ), which propagate according to Here the Green's functions G µν αβ (x, x , τ ekt , τ hydro ) describe the evolution and equilibration of perturbations from an early time τ ekt to a later time τ hydro .
Notably this linear response formalism provides a general framework to calculate the evolution of the energymomentum tensor, which requires limited microscopic input in the form of the nonequilibrium evolution of the background and linearized response functions. In this Letter, this microscopic input is calculated in QCD kinetic theory [10]. Starting from a microscopic gluon distribution function motivated by classical simulations of early time dynamics [17], we numerically solve the Boltzmann equation with the pure-glue leading order QCD collision kernels [10,30,31,33]. By analyzing the time dependence of the background distribution function and its perturbations, the evolution of the background energymomentum tensor and the Green's functions are then extracted from the kinetic theory simulations [33]. Equilibration time and conformal scaling.-Based on the approximate conformal symmetry of high-energy QCD, the rate of equilibration of the background energymomentum tensor T µν is governed solely by the QCD coupling constant λ = 4πα s N c along with a single dimensionful scale parameter. Since kinetic theory approaches hydrodynamics at late time, it is natural to express the dimensionful scale in terms of an asymptotic equilibrium quantity. Defining a pseudo-temperature T id ≡ (τ 1/3 T ) ∞ /τ 1/3 based on the asymptotic τ −1/3 dependence of the temperature in (conformal) Bjorken hydrodynamics, the equilibration rate is then determined by the kinetic relaxation time where the dependence on the coupling constant λ is encoded in the specific shear viscosity η/s. In contrast to earlier expectations [18] based on parametric estimates at (asymptotically) weak coupling, we find that for moderate values of η/s ( 1) the timescale τ hydro on which the nonequilibrium plasma approaches (viscous) hydrodynamics is determined by the equilibrium relaxation rate τ R [30,36]. Specifically, we observed in our kinetic theory simulations that the nonequilibrium evolution of the background energy density exhibits a universal scaling behavior. Expressing the evolution time τ in units of the relaxation time τ R , the ratio of the energy density e(τ ) to the asymptotic ideal hydrodynamic value e id. ∼ T 4 id , shown in Fig. 2, becomes independent of the microscopic coupling strength and smoothly interpo-lates between an approximate free-streaming behavior at early times and a hydrodynamic evolution at late times.
Based on the assumed symmetries, the universal curve in Fig. 2 determines the equilibration of the entire background energy-momentum tensor T µν (τ ). Such a collapse of the macroscopic evolution has also been referred to as a "hydrodynamic attractor" and is actively studied in the literature for different microscopic descriptions (mostly) in boost-invariant conformal systems [36][37][38][39]. Crucially, we found that the linear kinetic response functions G µν αβ (x, x , τ ekt , τ hydro ) used to propagate initial energy and momentum perturbations in Eq. (2) do not depend separately on the evolution time, background energy density, or the coupling constant, but only through the ratio of τ /τ R . One important consequence of this result is that the response functions need to be evaluated only once in a full kinetic theory simulation and can be reused for a different value of η/s or background temperature scale T id (τ ) by simple scaling transformations [33].
Besides its practical utility, the universality of the equilibration process in kinetic theory enables us to make pragmatic estimates of the time necessary before the plasma created in high-energy collisions can be described with hydrodynamics. We infer from Fig. 2 that the kinetic theory evolution approaches the hydrodynamic limit on time scales τ hydro ≈ 4πτ R . Relating the asymptotic constant (τ 1/3 T ) ∞ to the average entropy density per unit rapidity τ s , we obtain the following estimate for the hydrodynamization time

(3)
Here ν eff = ν eff (T eq ) is the effective number of degrees of freedom in the plasma at the equilibration temperature, as determined from the equilibrium relation s(T ) = ν eff 4π 2 90 T 3 . We use ν g = 2(N 2 c − 1) = 16 in our numerical simulations of a gluonic plasma, while for a realistic QCD equation of state, ν eff (0.4 GeV) ≈ 40 [40,41]. τ s for PbPb collisions √ s N N = 2.76 TeV is tightly constrained by experimental measurements and hydrodynamic simulations to be τ s ≈ 4.1 GeV 2 [31]. Dynamical description of pre-equilibrium stage:-We now present results obtained by applying our kinetic propagator KøMPøST to a realistic boost-invariant initial conditions of a central Pb-Pb collision at a centerof-mass energy √ s N N = 2.76 TeV. We start from an initial energy-momentum profile at τ ekt = 0.2 fm [42] given by the impact parameter dependent (IP)-glasma model [43,44], which provides a microscopic description of the classical Yang-Mills dynamics before the onset of equilibration. Subsequently, the locally averaged background energy-momentum tensor T µν for each point in the transverse plane is evolved in kinetic theory according to the universal evolution curve (see Fig. 2). Similarly, using appropriately scaled nonequilibrium response Average transverse and longitudinal pressures, PT = (T xx + T yy )/2 and PL = τ 2 T ηη , of a realistic heavy ion event evolved in succession by 2+1D Yang-Mills evolution (IP-glasma model) [43,44], QCD kinetic theory (KøMPøST) and relativistic viscous hydrodynamics [46][47][48].
In Fig. 3 we show overlapping theoretical descriptions of the evolution of the pressure anisotropy in the early stages of a realistic event. In the classical Yang-Mills field simulations the longitudinal pressure P L is initially negative, as is typical of a classical field configuration (cf. parallel plate capacitor with electric field E, where T ij = diag(E 2 , E 2 , −E 2 )/2 [50]). As the system evolves, the classical fields lose coherence and the longitudinal pressure approaches zero; the increasingly dilute system is then better described by kinetic theory. In the kinetic phase (τ ekt → τ hydro ) the energy-momentum tensor begins to equilibrate, such that ultimately the pressure approaches its equilibrium value of 1/3 of the energy density (for a locally equilibrated fluid of massless particles), up to corrections captured by viscous hydrodynamics. One clearly observes from Fig. 3 that the kinetic equilibration stage provides the missing link between the classical Yang-Mills evolution and the hydrodynamics, thus creating a self-consistent description of initial stages.
In order to illustrate the quality of the matching between the kinetic theory and the hydrodynamics, we investigate the robustness of the transverse energy and velocity profiles at τ out = 2.0 fm (in the hydrodynamic stage) under variations of the duration of the preequilibrium evolution τ hydro . In Fig. 4(a) we see that the energy density at τ out is essentially unchanged as the initialization time τ hydro is varied, indicating a smooth matching between the kinetic and hydrodynamic simulations. Figure 4(b) shows the corresponding plot for the transverse velocities. The transverse flow is also smoothly matched between the kinetic and hydrodynamic phases, with tension visible only at the edges of the fireball where our linearized approximation is pushed beyond its limits due to large gradients. We emphasize that KøMPøST also provides the viscous stress tensor π µν which is required for the subsequent (viscous) hydrodynamic evolution. We verified that the obtained π µν agrees reasonably well with the Navier-Stokes constitutive equations (π µν ≈ −ησ µν ), guaranteeing a consistent matching between the kinetic and hydrodynamic simulations [33].
Our current approach using kinetic response functions should be compared to preequilibrium modeling based on the long wavelength response [51][52][53] and free streaming [54,55]. We find that the first few terms in the long wavelength expansion capture most of the kinetic response [33]. However, to achieve a satisfactory level of agreement one needs to go beyond the leading order velocity response of Ref. [51]. Because of the universal-ity of the velocity response in conformal systems [31,51], we find that the velocity is well predicted even by free streaming. However, in the absence of longitudinal pressure, the free-streaming energy density decreases significantly slower than in kinetic theory or hydrodynamics. Moreover, the viscous stress tensor components π µν do not approach their hydrodynamic limit in the freestreaming description compromising the smooth matching to the hydrodynamic phase.
Conclusions & Outlook:-Nonequilibrium linear response theory captures the essential features of the early time preequilibrium dynamics of high-energy heavy ion collisions, and it provides a practical framework for connecting the theory of the initial state with the late time hydrodynamic evolution. Since the preequilibrium dynamics from kinetic theory matches smoothly onto viscous hydrodynamics, the approach naturally provides initial conditions for all second order hydrodynamic variables. There is no need to readjust the energy density, flow velocity or shear stress as the switching time τ hydro is changed, making the combined simulations increasingly predictive. Hadronic observables such as the multiplicity dN/dη or transverse momentum p T are now essentially insensitive to the switching time τ hydro [33].
Without the additional uncertainties from the initial state, near-thermal properties of the QCD medium can be increasingly constrained through multistage simulations and it would be interesting to reperform statistical inferences of transport coefficients η/s and ζ/s [56] with KøMPøST. Conversely, the properties of the earliest stage of the collisions can also be better inferred from experiments. For example, the color glass condensate predictions for the multiplicity as a function of centrality [57] should be revisited to take into account the rapid production of entropy during the kinetic pre-equilibrium phase.
A further benefit of handling the rich preequilibrium dynamics within kinetic theory is that the physics of chemical equilibration, early parton energy loss and electromagnetic emission can all be naturally described within the same QCD kinetic framework. The future inclusion of quark degrees of freedom in our initial state propagator KøMPøST [34] will be an important step in this direction.
Aside from improving the predictive power of stateof-the-art dynamical simulations of nucleus-nucleus collisions, the more complete understanding of the early time dynamics can also be used to assess the limits of applicability of the hydrodynamic description-a topic of considerable interest in light of recent experimental results in proton-proton (p+p) and proton-nucleus (p+A) collision [58]. Since a sufficiently long-lived hydrodynamic phase can be realized only when the equilibration time τ hydro is small compared to the system size R [59], one can use the estimate in Eq. (3) to determine the smallest system that can hydrodynamize. Since the entropy density per unit rapidity is approximately conserved during the later stages of the evolution, we can directly relate τ s at equilibrium time to the charged particle multiplicity dN ch /dη in the final state, according to τ s ≈ (S/N ch ) 1/A ⊥ dN ch /dη, where the transverse area of the system is estimated by A ⊥ ≈ πR 2 . Recalling that the ratio of entropy to charged particle number is S/N ch ≈ 7 for a hadron resonance gas [60,61], we obtain that the ratio τ hydro R 4π(η/s) 2 (4) is independent of the collision system size R and depends only on the charged particle multiplicity dN ch /dη [62] along with equilibrium and transport properties of the QGP. Then using Eq. (4) one concludes that in an optimistic scenario the minimum requirement for the hydrodynamic phase, i.e. τ hydro /R ≈ 1, is reached if the charged particle multiplicity is at least dN ch /dη 8 for a small η/s = 1/(4π). However, this estimate is rather sensitive to the transport properties of the QGP as for a larger value of specific shear viscosity η/s = 2/(4π) the minimal required multiplicity increases to dN ch /dη 63. Based on the experimentally measured minimum bias multiplicities of dN ch /dη 6 in 7 TeV p+p [63] and dN ch /dη 16 in 5.02 TeV p+Pb collisions [64], one concludes that a hydrodynamic phase is unlikely to emerge in most p+p collisions. However in p+Pb collisions, where events with several (up to ∼ 10) times the minimum bias multiplicity can be observed, it seems possible that a hydrodynamically flowing QGP can be formed.
The authors would like to thank Björn Schenke for the 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, and 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