Initial conditions for hydrodynamics from weakly coupled pre-equilibrium evolution

# Initial conditions for hydrodynamics from weakly coupled pre-equilibrium evolution

Liam Keegan Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway    Aleksi Kurkela Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway    Aleksas Mazeliauskas Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway    Derek Teaney Faculty of Science and Technology, University of Stavanger, 4036 Stavanger, Norway
###### Abstract

We use effective kinetic theory, accurate at weak coupling, to simulate the pre-equilibrium evolution of transverse energy and flow perturbations in heavy-ion collisions. We provide a Green function which propagates the initial perturbations to the energy-momentum tensor at a time when hydrodynamics becomes applicable. With this map, the complete pre-thermal evolution from saturated nuclei to hydrodynamics can be modelled in a perturbatively controlled way.

###### Keywords:
Quark-Gluon Plasma, Perturbative QCD, Heavy Ion Phenomenology
\preprint

CERN-TH-2016-117 Theoretical Physics Department, CERN, Geneva, Switzerland Department of Physics and Astronomy, Stony Brook University, Stony Brook, New York 11794, USA \arxivnumber1605.04287

## 1 Introduction

Viscous relativistic hydrodynamics provides a remarkably detailed and phenomenologically successful description of the expansion of the Quark Gluon Plasma (QGP) in the ultrarelativistic heavy-ion collisions realized at the BNL Relativistic Heavy Ion Collider (RHIC) and at the CERN Large Hadron Collider (LHC) Heinz:2013th (); Luzum:2013yya (); Teaney:2009qa (). Hydrodynamics is an effective theory based on an assumption that the medium is sufficiently close to local thermal equilibrium that the full stress tensor can be expanded in gradients of the energy and momentum densities Baier:2007ix (). However, due to the singular geometry of heavy ion collisions, the gradients diverge at early times, and the hydrodynamic approach does not apply during the initial stages of the evolution. Indeed, hydrodynamic simulations start at some sufficiently late initialization time , when the gradient expansion becomes a useful approximation scheme. The initial conditions for hydrodynamics at are generally unknown, and must be parametrized and fitted to data Bernhard:2015hxa (). This procedure often neglects any prethermal evolution, and limits the empirical determination of the transport coefficients of the QGP Liu:2015nwa ().

A useful prethermal model should smoothly and automatically approach hydrodynamics. If this is the case, the combined pre-thermal and hydrodynamic evolutions will be independent of the initialization time vanderSchee:2013pia (); Romatschke:2015gxa (); Kurkela:2016vts (). In most simulations the prethermal evolution is either completely neglected Niemi:2015qia (), or modelled in a way that does not contain the correct physics to produce hydrodynamic flow Broniowski:2008qk (); Schenke:2012wb (); Liu:2015nwa (). In addition, in some models (such as the successful IP-glasma model Schenke:2012wb () motivated by parton saturation) the initial conditions contain strong gradients which limit the effectiveness of the hydrodynamic derivative expansion Niemi:2014wta (); Noronha-Hostler:2015coa (). Different hydrodynamic codes regulate these extreme initial conditions in different ad hoc ways, e.g. by arbitrarily setting the shear stress tensor to zero when the hydrodynamics is initialized. Again, these ambiguities limit the ability of hydrodynamic simulations to determine the transport properties of the QGP.

In the limit of weak coupling the approach to hydrodynamics, or hydrodynamization, is described by an effective kinetic theory (EKT) Arnold:2002zm (), which takes into account the non-trivial in-medium dynamics of screening and the Landau-Pomeranchuk-Migdal suppression of collinear radiation. In ref. Kurkela:2015qoa () (which includes one of the authors), it was shown that the EKT, starting with initial conditions motivated by the Color-Glass Condensate (CGC) saturation framework Iancu:2002xk (); Iancu:2003xm (); Gelis:2010nm (); Gelis:2007kn (); Lappi:2011ju (), reaches hydrodynamics in a phenomenologically reasonable time scale of , where is the (adjoint representation) saturation scale, which is estimated to be of order of few GeV for central heavy-ion collisions at the LHC. This first calculation used the EKT to monitor the equilibration of a uniform plasma of infinite transverse extent during a Bjorken expansion.

Transverse gradients in the profile will initiate flow during the equilibration process. This preflow and the accompanying modifications of the initial energy density profile will influence the subsequent hydrodynamic evolution. The goal of the current paper is to use the EKT to precisely determine the preflow and the components of the energy momentum tensor that should be used to initiate the hydrodynamic evolution for a given initial energy density profile. Although the kinetic theory calculation can be used to match different models for the initial energy profile to hydrodynamics, the weak coupling approximations made in the IP-glasma model lead naturally to effective kinetic theory.

Figure 1 shows a typical transverse (entropy) profile that is used in

current hydrodynamic simulations Mazeliauskas:2015vea (). Clearly during the equilibration process the profile will change and generate initial flow. The equilibration time, , is short compared to the nuclear radius, . For this reason the prethermal evolution is insensitive to the global collision geometry. Indeed, we may decompose the transverse plane into causally disconnected patches of size whose prethermal evolution can be separately determined. In these patches, the global nuclear geometry determines a small gradient that can be considered as a linear perturbation over a translationally invariant background. Thus, corrections to initial conditions for hydrodynamics from the global geometry are of order  Vredevoogd:2008id (). In addition to the global geometry, the initial energy density profile includes event-by-event fluctuations at smaller scales set by the nucleon size , which is comparable to the causal horizon . Event-by-event fluctuations at these length scales are suppressed by where is the number of participating nucleons in the event, . Therefore, such fluctuations can also be treated in a linearized way as fluctuations over a translationally invariant background. The structure of the initial profile at even smaller scales is less well known, but in models based on CGC, one expects fluctuations to subnuclear scales of order the saturation momentum, .

Finally, an important scale is set by the mean free path, which in a weakly coupled theory is of order for states not too far from equilibrium. In practice, this length scale is comparable, though slightly shorter than the causal horizon and the nucleon scales. Without the scale separation, the medium prethermal response to initial perturbations in the transverse plane can only be computed by a calculation within the EKT. Fortunately, as discussed above linearized kinetic theory is sufficient to determine this response.

To summarize, our strategy is to use linearized kinetic theory to follow the hydrodynamization of perturbations on top of a far-from-equilibrium Bjorken background with translational symmetry in the transverse directions. This determines the stress tensor for hydrodynamics at the initialization time. The length scales of relevance are the nuclear-geometry, the nucleonic scale, the causal horizon , and the mean free path

 R≫Rp∼cτinit∼1λ2Teff. (1)

By linearizing the problem and solving for the response, we will determine a Green function describing how an energy fluctuation at the earliest moments, , evolves during the equilibration process to the hydrodynamic fields, i.e. the energy and momentum densities, and respectively. We will verify that the subsequent evolution is described by second order hydrodynamics to certifiable precision.

In section 2 we outline the linearized EKT, and study the linear response of the EKT in equilibrium. In section 3 we systematically study the approach to equilibrium of Fourier modes of specified , starting with a far from equilibrium initial state. In section 4 we Fourier transform these results and determine a coordinate space Green function which produces the appropriate initial conditions for hydrodynamics at when convolved with a specified initial state. We also analyze the long wavelength limit of these Green functions, making contact and providing additional insight into previous work on preflow Vredevoogd:2008id (). Finally, we discuss our conclusions in section 5.

## 2 Linearized kinetic theory

### 2.1 Setup

At weak coupling the non-equilibrium evolution of the boost invariant color and spin averaged gluon distribution function is described in terms of an effective kinetic equation Arnold:2002zm ()

 ∂τfx⊥,p+p|p|⋅∇x⊥fx⊥,p−pzτ∂pzfx⊥,p=−C[fx⊥,p], (2)

where the effective collision kernel incorporates the elastic and inelastic processes as required for a leading order description in the coupling constant , which is the only parameter of the EKT. The kinetic theory is valid when the occupancies are perturbative and when the relevant distance scales are larger than the typical Compton wavelength of the particles . The details of the scattering kernel have been discussed in refs. Arnold:2002zm (); Kurkela:2015qoa (); Keegan:2015avk () and are briefly repeated here in the appendix A. We use the isotropic screening approximation from Kurkela:2015qoa () which is leading order accurate for parametrically isotropic systems . (Here and below and denote the longitudinal and transverse pressures.) In the current paper we will consider only gluonic degrees of freedom and assume that the contribution of quarks is suppressed during the pre-equilibrium evolution.111The initial far-from-equilibrium state is parametrically dominated by gluons. Once the plasma has thermalized it should contain also fermionic degrees of freedom. However, the production of fermions is suppressed by larger color factors , and by Pauli blocking factors (while scattering of gluons is Bose enchanced). It is therefore plausible that the system hydrodynamizes before it is chemically equilibrated.

We split the distribution function into a translationally symmetric background and a linearized perturbation with a wavenumber in the transverse plane

 fx⊥,p=¯fp+∫d2k⊥(2π)2A(k⊥)eik⊥⋅x⊥δfk⊥,p, (3)

where characterizes the initial density profile. The kinetic equations for the background and the (complex) fluctuation then read

 (∂τ−pzτ∂pz)¯fp = −C[¯f], (4a) (∂τ−pzτ∂pz+ip⊥⋅k⊥p)δfk⊥,p = −C[¯f,δf], (4b)

where is the collision kernel linearized in (see appendix A for details).

### 2.2 Hydrodynamization close to equilibrium

Before studying the equilibration process, we will analyze the linear response of the EKT close to equilibrium, corresponding to the limit of (4). Our goal in this section is to determine at what wavenumbers (characterized by ) linearized energy-momentum perturbations are described by hydrodynamics for an equilibrated background.

The dispersion relation for the sound mode to second order in the hydrodynamic expansion reads Baier:2007ix ()

 ω=csk−i43ηe+pk2+43ηe+p(csτπ−23csηe+p)k3, (5)

where for conformal equation of state and , are known transport coefficients at weak coupling Arnold:2003zc (); York:2008rr (). For (corresponding to ) the hydrodynamic coefficients read , , and . We will quantify at what numerical values of the corrections to (5) become sizeable.

To this end, the kinetic theory is initiated in local thermal equilibrium with a spatially varying temperature, , and corresponding phase space distribution

 δf(1)k⊥,p=−δTTp∂p¯fp, and¯fp=1ep/T−1. (6)

Gradients in the energy density drive momentum perturbations in due course, and the frequency of the subsequent (damped) oscillations determines the real part of the dispersion relation. Numerically, we obtained the oscillation frequency by measuring the time interval between the successive nodes.

The results for various are depicted in figure 2 for . We see that at small , , the dispersion relation is well described by the result of ideal hydrodynamics, and the approach to ideal hydrodynamics is described by the 2nd order corrections of (5) (note that the real part of the frequency does not get a first order correction). Indeed, for , the second order hydrodynamic theory matches well with the EKT.

At higher values of , the EKT finally saturates at in contrast to the strict (unresummed) second order hydrodynamics. For this happens only at rather large values of , . For these wavenumbers, the wavelength of the perturbation is comparable to the typical gluon Compton wavelength, and the linear response of the system cannot be reliably computed with kinetic theory in this regime.

We conclude that for , the smallest scales that hydrodynamize have . Varying the value of (not shown), we find that the scale where hydrodynamics breaks downs tracks the shear viscosity, with varying . Note that for smaller , the saturation to can take place within the regime of validity of the effective theory.

## 3 Hydrodynamization of fluctuations far from equilibrium

We now move on to study the hydrodynamization of spatially dependent fluctuations on top of a far-from-equilibrium boost invariant background. As discussed in Kovner:1995ja (); Lappi:2006fp (); Gelis:2013rba (); Kurkela:2015qoa (), at very early times the far-from-equilibrium gluonic system in the midrapidity region is parametrically over-occupied , and the dynamics is described with coherent classical gauge fields rather than with particles. This part of the evolution is characterized by negative values of the longitudinal pressure , which is a result of the coherence of the approximately boost invariant fields. However, classical numerical simulations Lappi:2011ju (); Gelis:2013rba () (as well as analytical series solutions to the classical equations of motion Chen:2015wia (); Li:2016eqr ()) show that in a timescale , the coherence is lost, the longitudinal pressure approaches zero , and the occupancies become perturbative Kurkela:2011ti (); Kurkela:2011ub (); Berges:2013fga (); Baier:2000sb (). At this point the system may be passed to the EKT Mueller:2002gd (); Jeon:2004dh (); York:2014wja (); Kurkela:2015qoa ().

Following Kurkela:2015qoa (), we take as our initial condition at a parametrization

 f(pz,p⊥) =2λAf0(pzξ/p0,p⊥/p0), (7) f0(^pz,^p⊥) =1√^p2⊥+^p2ze−2(^p2⊥+^p2z)/3, (8)

where , , and . The parameters and are motivated by classical simulations where and . The amplitude, , is adjusted so that energy per rapidity

 τ0e(τ0)=τ0νg∫d3p(2π)3|p|f(pz,p⊥), (9)

matches the results of classical simulations Lappi:2011ju (), where

 τ0e(τ0)≃0.358τ0νgQ4sλ. (10)

Here is the number of gluonic degrees of freedom. With these parameters, the number of gluons and the mean in the EKT at are

 dNd2x⊥dy= 0.232νgQ2sλ,√⟨p2T⟩=1.8Qs, (11)

which roughly matches the classical Yang-Mills simulations.

We will follow the response to two specific initial perturbations of these initial conditions, which provide an independent basis for describing energy and momentum fluctuations in the transverse plane. For the energy fluctuation, we take

 δf(1)k⊥,p(τ0) =−δQsQsp∂p¯fp, (12)

which results from varying the saturation scale in (7), . This initial condition is studied in the bulk of the paper. For the momentum fluctuation we take

 δf(2)k⊥,p(τ0) =−δQsQs^k⊥⋅∂p¯fp, (13)

which does not perturb the local energy density. The response to this initial condition is recorded in appendix C.

Without loss of generality, we can choose the wave vector to point in -direction. Then at any time, the energy and momentum perturbations are defined as

 δe(τ,k) ≡δT00=νg∫d3p(2π)3p0δf, (14) gx(τ,k) ≡δT0x=νg∫d3p(2π)3pxδf, (15)

and their evolution is governed by the linearized conservation equations

 ∂τe(τ)=−e(τ)+Tzz(τ)τ, (16a) ∂τδe(τ,k)+ikgx(τ,k)=−δe(τ,k)+δTzz(τ,k)τ, (16b) ∂τgx(τ,k)+ikδTxx(τ,k)=−gx(τ,k)τ, (16c)

where is the energy-momentum tensor perturbation caused by . If the system is described by hydrodynamics then , and are determined through the constitutive equations by the first moments of the particle distribution function , and . For conformal second order viscous hydrodynamics these relations are

 Tzz(τ) =13e−43ητ−89τπη−λ1τ2, (17a) δTxx(τ,k)= δe(τ,k)e[13e+13ητπk2+12τη−2(λ1−ητπ)9τ2] −ikgx(τ,k)e[η−1τ(η22e+ητπ2−23λ1)] (17b) δTzz(τ,k)= δe(τ,k)e[13e−16ητπk2−1τη+4(λ1−ητπ)9τ2] +ikgx(τ,k)e[12η−1τ(η24e+23λ1)], (17c)

where the constitutive equations for ideal (or first order viscous) hydrodynamics can be recovered by setting (or ).

### 3.1 Evolution of the background energy density

Before studying the perturbations, we will study the equilibration of the background energy density, elaborating on the original study Kurkela:2015qoa (). In figure 3(a) we compare the energy momentum tensor combination in the kinetic theory simulation to the constitutive equation, (17a). The combination is motivated by the conservation law in (16a).

At early times the system evolves approximately according to free streaming222During this part of the evolution, the system evolves according to the nonthermal attractors discussed in e.g. Baier:2000sb (); Kurkela:2011ub (); Berges:2013fga (); Kurkela:2015qoa (). The nonthermal attractors are characterized by , and for the current discussion the fine details of the attractor are irrelevant, and the evolution resembles that of free streaming., with and . As already noticed in Kurkela:2015qoa (), the constitutive equations give an increasingly accurate description of the EKT stress tensor as a function of time. While the ideal constitutive equations are rather far from the EKT at all relevant times, the viscous and 2nd order equations quickly converge to the EKT. Note that an accidental (approximate) cancellation of makes the second order correction anomalously small Luzum:2008cw (), and only at rather late times after non-hydrodynamic modes have almost completely decayed does second order hydrodynamics finally improve the first order result (not shown). By times , the viscous constitutive equations agree with the EKT within .

It is noteworthy that the EKT interpolates smoothly between the free streaming and viscous hydrodynamic evolutions without an extended period during which the evolution is not approximately described by one or the other approximation scheme. At the evolution is somewhat closer to free streaming, and using hydrodynamics at this point is a rough, though perhaps acceptable, approximation. At hydrodynamics is a better approximation, but at this time the causal horizon is becoming comparable to the nuclear radius.

Given the agreement with the constitutive equations, one can use the hydrodynamic equations to propagate the system forward in time. At late times ideal hydrodynamics is valid, and the entropy per area per rapidity approaches a constant

 limτ→∞τs(τ)≡νgΛ2sλ. (18)

(Here the parametrization is motivated by the scaling of the initial multiplicity with the saturation scale in (10).) Dimensional reasoning indicates that is proportional to . Taking the data presented in figure 3(a) we may extrapolate to determine the proportionality coefficient

 Λ2s=1.95Q2s. (19)

Here we have used the ideal equation of state to convert energy density to entropy density. Since the entropy per gluon of an ideal gluon gas is , (19) implies the asymptotic number of gluons per area per rapidity is more than a factor of two larger than the input number of gluons at

 dNd2x⊥dη∣∣∣τ→∞=2.33dNd2x⊥dη∣∣∣τ=τ0. (20)

At the entropy and gluon multiplicity have reached only {72, 82}% of their asymptotic values, corresponding to gluon multiplication factors of respectively.

Finally, let us make phenomenological contact with more complete hydrodynamic simulations of heavy ion collisions, and estimate the saturation momentum required by phenomenology. The initial entropy in hydrodynamics is normally adjusted to reproduce the mean multiplicity. Using the computer code from one such hydrodynamic simulation at the LHC Mazeliauskas:2015vea (), we computed the average entropy per area at the hydrodynamic initialization time333Specifically, for the event-by-event hydro code described in ref. Mazeliauskas:2015vea () we first created a smooth entropy density profile, , by averaging over events for a 0-5% centrality class, . This average is shown in figure 1. We then computed a single averaged entropy density by averaging with as a radial weight.

 ⟨τinits(τinit)⟩=4.13GeV2. (21)

Entropy production during the subsequent hydrodynamic evolution is small, approximately 15%, and therefore this constant is approximately independent of the initialization time. With (18) and (19), setting in the EKT will roughly reproduce the entropy in hydrodynamic simulations provided the system is passed to hydrodynamics at where the entropy in the EKT has reached 70% of its asymptotic value in (18).

In figure 3(b) we present the time evolution of the effective temperature (determined from the energy density and the equation of state) in physical units for . Its time dependence at late times is well described by first order viscous hydrodynamics with the asymptotic value

 limτ→∞(T+23ηsτ)τ1/3=0.763Q2/3s. (22)

The temperature at is , or close to three times the pseudo-critical temperature. At these temperatures, modern weakly coupled techniques can be expected to work reasonably, justifying our approximation scheme. Similarly, with this value of an initialization time of corresponds to . The elliptic and triangular flows in central events develop on a later time scale, of order , and therefore an initialization time of this order may be acceptable. A more complete study including fermions will be needed to definitively answer this question.

### 3.2 Evolution of the perturbations

We now move on to describe the evolution of the linearized perturbations on top of the thermalizing non-equilibrium background. We first follow the evolution of an energy density perturbation, which we start with the initial condition of (12) with different values of . In figures 4(a)-(d) we show -component of perturbation energy momentum tensor compared with ideal, viscous, and second order hydrodynamic constitutive equations of (17b). The lines have been normalized by the background energy density , so that any observed damping is due to nontrivial dynamics associated to the spatial inhomogeneity.

We consider fixed values of , which do not correspond to fixed values of , as the effective temperature is changing due to the expansion (see figure 4). For very large wavelengths with , we observe that while ideal constitutive equations have rather large corrections, these are well accounted for by the viscous and 2nd order equations, roughly at the same time as the background constitutive equation is satisfied, i.e. at times after . Determining the temperature of the background from Landau matching condition (see figure 3 (right)), we find that at times , the wavelength in units of temperature is . As discussed in section 2.2 and in figure 2 these values of are accurately described by second order hydrodynamics.

We see that even for larger the hydrodynamic constitutive relations are approximately fulfilled when the background has hydrodynamized around . However, larger values of correspond to larger values of , and even at late times there are corrections to the constitutive equations. While these corrections are moderate for for which , they remain for . It is therefore questionable whether it is justifiable to pass these short scales to hydrodynamic description at any time.

## 4 A Green function for hydrodynamics

We now move on to describe how the response to the linearized perturbations in EKT can be used in a hydrodynamic simulations to encapsulate the far-from-equilibrium dynamics of transverse perturbations during the time scales between and .

In order to construct the initial state for hydrodynamics at from a given geometry at , the linear response of the components of to the initial perturbation are needed. The constitutive relations reduce the number of independent components of the energy momentum tensor, so it suffices to specify only and . Figure 5(a) displays the energy and momentum response functions ( and respectively) to an initial energy perturbation in -space

 δe(τ,k)e(τ) ≡~E(k;τ,τ0)δe(τ0,k)e(τ0), (23) gx(τ,k)e(τ) ≡−i~G(k;τ,τ0)δe(τ0,k)e(τ0). (24)

The results are presented at two suggested initial times . We will analyze the response functions at asymptotically small and in coordinate space in the next two subsections. Analogous results for the response to an initial momentum perturbation are given in the appendix C.

### 4.1 The kinetic theory response at asymptotically small k

The most important contribution to the flow arises from the average nuclear geometry, which is smooth on spatial scales of order . For this reason the flow due to the average geometry is determined by the limit of the response functions. This section will provide an analytic understanding of this limit, i.e. the intercept of and the slope of in figure 5.

First, we will determine how the long wavelength energy perturbations in the transverse plane change as a function of time. Returning to the conservation equations, (16a) and (16b), and setting , we have

 ∂τe(τ)=−e(τ)+Tzz(τ)τ, (25) ∂τδe(τ)=−δe(τ)+δTzz(τ)τ. (26)

From these equations the fractional perturbations in the transverse plane remain constant in time in the free streaming limit (where and are zero), and in the hydrodynamic limit (where and are one third and ). Outside of these limits is not constant in time.

However, a constant of the motion at can be constructed whenever the hydrodynamic gradient expansion is applicable. Indeed, by dimensional analysis, an all order constitutive equation at must take the following form

 Tzz=ef(e1/4τ), (27)

where is an order one function and . Then straightforward steps show that to all orders in the gradient expansion

 limk→0δe(τ,k)3e−Tzz=const. (28)

For conformally invariant theories with this can be written as

 limk→0δe(τ,k)e+Txx=const. (29)

In figure 6(a) we present the time evolution of and relative to their initial values. For our initial conditions remains very nearly constant throughout the entire evolution. Using this result, the change in can be determined by the ratio of at the initial and final times, when is approximately and respectively. This reasoning leads to an asymptotic relation between the initial and final energy perturbations

 limτ→∞δe(τ)e(τ)=89δe(τ0)e(τ0), (30)

which is shown in figure 6(a).

Next, we will determine the velocity at (asymptotically) small as a function of time from the pre-thermal evolution. From the conservation equations for perturbations, (16b) and (16c), the momentum perturbations at small satisfy

 ∂τ(τgxik+12δeτ2)=−τ2(2δTxx+δTzz−δT00). (31)

For conformal theories with the right hand side of (31) is zero and

 τgxik+12δeτ2=const. (32)

At late times and in coordinate space this condition reads

 T0x(τ)T00(τ)=−12τ∂xT00(τ)T00(τ), (33)

which was first noted in Vredevoogd:2008id (). Here we have shown that this relation is a consequence of conformal symmetry (see also Vredevoogd:2008id ()) and the small limit.

Using (32) and the definition , the velocity as a function of time is given by

 vxik=−τ2δee+Txx(1−δe(τ0)τ20δe(τ)τ2). (34)

Thus, after a brief transient period of order , the velocity is directly proportional to time

 vx=τ2(−∂xe(τ,x)e(τ)+Txx(τ)),−∂xe(τ,x)e(τ)+Txx(τ)=const. (35)

In figure 6(b) we compare the growth of the velocity with time given by (34) with a simple estimate based on (35). The simple estimate does a remarkably good job for all times.

### 4.2 Response in coordinate space

To construct the initial conditions for hydrodynamics with the correct prethermal evolution, we determine the Green functions and which convert the initial profile of energy perturbations to the required energy and momentum fluctuations at thermalization time

 δe(τ,x)e(τ) =∫d2x′δe(τ0,x′)e(τ0)E(|x−x′|;τ,τ0), (36a) gi(τ,x)e(τ) =∫d2x′δe(τ0,x′)e(τ0)(x−x′)i|x−x′|G(|x−x′|;τ,τ0). (36b)

Currently hydrodynamic simulations often smooth the initial conditions before starting the hydrodynamic evolution by convolving the energy density with a Gaussian444See ref. Noronha-Hostler:2015coa () for a current discussion of the observables that are influenced by this arbitrary regulator.. (36) by contrast smooths the initial conditions in a physical way, and provides an attractive alternative to this ad hoc procedure.

The EKT is applicable for distance scales that are larger than the Compton wavelength of the particles . This limits the accuracy of the Green function in spatial domain that can be reached in a computation based on kinetic theory. In order to fold this uncertainty into our result, we regulate our Green function by convoluting with a Gaussian weight, , with , and with a width of the order of the initial Compton wavelength . In momentum space this corresponds to suppressing the large contributions by an exponential envelope

 E(|x|;τ,τ0) =∫d2k(2π)2eik⋅xe−σ2k2/2E(|k|;τ,τ0), (37) G(|x|;τ,τ0) =∫d2k(2π)2(−i^k⋅^x)eik⋅xe−σ2k2/2G(|k|;τ,τ0). (38)

The regulated Green functions are shown in figure 7 at the initialization times (for details of the Fourier transform see appendix B.) At the system has spent a significant proportion of the total evolution time with small longitudinal pressure , and therefore the resulting response is similar to the free streaming prediction (see appendix B). However, the Green function in figure 7(a) is peaked for , suggesting a slight deflection from the free streaming trajectory. Additionally, the energy perturbation is negative at small , which is indicative of a nascent approach to hydrodynamics. At later times, such as in figure 7(b), these differences become more pronounced. Similar features are visible in the momentum response to an initial energy perturbation shown figure 7(c) and (d). Finally, in figure 8 we show the Green functions at later times and , and compare to linearized second order hydrodynamics ((16) and (17)) with initial conditions taken from the results. Between and , the hydrodynamics overdamps the high modes (see also Romatschke:2015gic ()), and the response is broader than the predictions of kinetic theory. However, these Green functions will be convolved with the initial conditions, and thus the resulting hydrodynamic initial state is mostly sensitive to the first moments of these kernels. The moments of the EKT and hydro kernels are determined by the small behaviour of the response functions, which agree to a few percent (not shown). At later times , the response is largely determine by Fourier modes in the hydrodynamic regime , and the EKT and hydro kernels are visually similar.

To summarize, the hydrodynamic evolution sets in early at rather large anisotropies, and the hydrodynamic constitutive equations are approximately satisfied as soon as the starts to significantly deviate from the free streaming expectation, . For this reason the time interval when the evolution is not described by free streaming or hydrodynamics is comparatively brief (see figure 3(a)), and as hydrodynamics becomes marginally applicable at , the Green function closely resembles the free streaming result. Therefore, an approach where the evolution is described by free streaming until seems well motivated Broniowski:2008qk (); Liu:2015nwa (), provided that the correct value of is used. However, such an ad hoc approach does not account for some of the qualitative details of the Green function, such as the depletion of the energy density in the interior region.

## 5 Discussion

In this paper, we have provided a bridge between the far-from-equilibrium initial conditions of heavy-ion collisions and hydrodynamized plasma. Our main result is the coordinate space Green functions (see (36) and figure 7), which can be used to filter the pre-equilibrium energy density to find the full energy-momentum tensor for hydrodynamics at the initialization time. The procedure can be implemented in complete hydrodynamic simulations, removing one source of uncertainty. Perhaps more importantly, the approximations in the EKT are compatible with the IP-Glasma setup, and thus the whole evolution from saturated nuclei to hydrodynamics can be comprehensively modelled within a perturbatively controlled framework.

We provide the coordinate space Green functions at two different suggested initialization times, . At the earlier initialization time, , there are significant (though bearable) corrections to the constitutive relations due to non-hydrodynamic modes (see figure 3(a) and figure 4). By the constitutive relations at small are well satisfied, and the subsequent evolution is reasonably captured by second order hydrodynamics555When examining figure 8, one must remember that the full Green functions will be convolved with the initial conditions, and thus the response of the system is mostly sensitive to the first moments of the kernels in figure 8. The first EKT moments (i.e. the small behavior of the response) agree with the hydrodynamics to the percent level. (see figure 8). The approximate overlap of the two 2nd order viscous lines in figure 8, which correspond to initializing the hydro at , demonstrates that the subsequent hydrodynamical evolution is indeed rather insensitive to the initialization time. In section 4 we examined the (asymptotically) small limit of the Green functions, and confirmed (and clarified) a preflow estimate by Vredevoogd and Pratt Vredevoogd:2008id () (see figure 6).

The hydrodynamics that the EKT follows is characterized by the weak coupling value of  Arnold:2003zc (), which is significantly higher than the AdS/CFT result  Policastro:2001yc (), and current phenomenological estimates, which assume that is independent of the temperature. Recent analyses have relaxed the temperature independence of , and shown that the value of at higher temperatures is poorly constrained by data Niemi:2015qia (). Since it is the high temperature regime that is most relevant for the transition to hydrodynamics, we believe that the current kinetic theory results for the initial stages can be consistent with hydro phenomenology, provided decreases towards the strong coupling result as the system cools towards .

Nevertheless, to apply our results to a hydro simulation with lower viscosity than the perturbative expectation, we note that the two initialization times correspond in units of the hydro parameters to . The scaled times can be used to initialize simulations when the transport coefficients differ. Such an approach is supported by the reasonably good scaling properties of the hydrodynamization times and prethermal evolution as function of the coupling constant when expressed in terms of the hydrodynamic variables Kurkela:2015qoa (); Keegan:2015avk ().

Although our EKT description can be further improved by inclusion of fermionic degrees of freedom and by improving the connection to the early classical evolution, we believe that it already provides a physically sound picture of the approach to hydrodynamics and can be used to initialize all components of the energy-momentum tensor for subsequent hydrodynamic evolution. This eliminates a source of uncertainty in current simulations, and provides a satisfyingly complete description of the early time evolution in heavy ion collisions.

\acknowledgments

The authors thank François Gelis, Keijo Kajantie, Tuomas Lappi, Risto Paatelainen, Jean-François Paquet, and Urs Wiedemann for useful discussions. The simulations were performed using computer resources of UiS and high-performance LIred computing system at the Institute for Advanced Computational Science at Stony Brook University. A.M. and D.T. work was supported in part by the U.S. Department of Energy under Contracts No. DE-FG-88ER40388. Finally, A.M. and D.T would like to thank CERN Theoretical Physics Department for the hospitality during the short-term visit.

## Appendix A Collision kernel

In this appendix we provide additional details on the collision kernels used in (4). The collision kernel for the uniform background contains terms arising from elastic scatterings and inelastic collinear splittings

 C[f]=C2↔2[f]+C1↔2[f]. (39)

The two collision terms read  Arnold:2002zm (); Kurkela:2015qoa (); Keegan:2015avk ()

 C2↔2[f](p) =14|p|νg∫d3k2k(2π)3d3p′2p′(2π)3d3k′2k′(2π)3|M(p,k;p′,k′)|2(2π)4δ(4)(P+K−P′−K′) ×{fpfk[1+fp′][1+fk′]−fp′fk′[1+fp][1+fk]} (40)

and

 C1↔2[f](p) =(2π)32|p|2νg∫∞0dp′dk′δ(|p|−p′−k′)γ(p;p′^p,k′^p)×{fp[1+fp′^p][1+fk′^p]−fp′^pfk′^p[1+fp]} (41)

where is the unit vector parallel to , and capital letters denote null 4-vectors, i.e. . The effective elastic and inelastic scattering matrix elements contain non-trivial structures arising from the soft and collinear divergences, which are dynamically regulated by the in-medium physics.

For the most of kinematics the effective elastic scattering element is given by666 Equations (42) and (46) have some minor typos corrected compared to refs. Kurkela:2015qoa (); Keegan:2015avk ().

 |M|2=2λ2νg(9+(s−t)2u2+(u−s)2t2+(t−u)2s2). (42)

For a soft gluon exchange with the momentum transfer in -channel (or in -channel) the collision matrix is proportional to , and thus suffers from a soft Coulomb divergence. It is regulated by replacing

 q2t→(q2+2ξ20m2)t, (43)

in the denominators of divergent terms (similarly for the -channel). Here is the thermal asymptotic mass of the gluon defined as

 m2=2λ∫d3p(2π)3fp|p|. (44)

The coefficient is fixed so that the matrix element reproduces the drag and momentum diffusion properties of soft scattering at leading order for isotropic distributions York:2014wja ().

 γ(p^p;p′^p,k′^p)=p4+p′4+k′4p3p′3k′3νgλ8(2π)4∫d2h(2π)22h⋅ReF, (45)

where the equation for accounts for splitting due to multiple scatterings with transverse momentum exchange , and momentum non-collinearity

 2h= iδE(h)F(h)+λT∗2∫d2q⊥(2π)2A(q⊥) (46) ×[3F(h)−F(h−p′q⊥)−F(h−k′q⊥)−F(h+pq⊥)].

with and . In the isotropic screening approximation

 A(q⊥)=(1q2⊥−1q2⊥+2m2). (47)

Both and are self-consistently evaluated at each time step.

The linearized collision kernels are obtained trivially by replacing in the integrands of A and A and linearizing in . In addition one has to take into account the linear variation of the thermal mass and the effective temperature in the scattering matrix elements (42) and (45)

 δm2 =2λ∫d3p(2π)3δfp|p|, (48) δT∗ =λm2∫d3p(2π)3δfp(1+2fp)−δm2m2T∗. (49)

where and are evaluated from the unperturbed background distribution.

## Appendix B Fourier transform of Green functions

Here we provide details of performing Fourier transforms in (37) and (38) to obtain spatial Green functions shown in figures 7 and 8. The two dimensional Fourier transforms can be straightforwardly reduced to one dimensional Hankel transforms

 E(|x|;τ,τ0) =∫∞0dk2πk~E(k;τ,τ0)e−σ2k2/2J0(k|x|), (50) G(|x|;τ,τ0) =∫∞0dk2πk~G(k;τ,τ0)e−σ2k2/2J1(k|x|). (51)

Integrals in (50) and (51) were done numerically by using cubic interpolation for and within the available range of wavenumbers To avoid the oscillatory behaviour due to a sharp cut-off at , we extrapolated the Green functions until the Gaussian envelope smoothly cuts off the integral. For extrapolation at large we used functional forms motivated by free streaming results: and , where coefficients