Isotropization and hydrodynamization in weakly coupled heavy-ion collisions

We numerically solve 2+1D effective kinetic theory of weak coupling QCD under longitudinal expansion relevant for early stages of heavy-ion collisions. We find agreement with viscous hydrodynamics and classical Yang-Mills simulations in the regimes where they are applicable. By choosing initial conditions that are motivated by color-glass-condensate framework we find that for Q=2GeV and $\alpha_s$=0.3 the system is approximately described by viscous hydrodynamics well before $\tau \lesssim 1.0$ fm/c.


INTRODUCTION
In the weak coupling picture of the pre-thermal evolution of heavy-ion collisions, the post collision debris that end up in the mid-rapidity region undergo several stages that are characterised by widely different physics and degrees of freedom. On the one hand, according to the saturation paradigm [1], at very early times τ ∼ Q −1 s , where Q s is the typical energy scale right after the collision, the energy is deposited in strong color fields. These strong fields admit a description in terms of classical Yang-Mills theory to leading order in 't Hooft coupling λ = 4πN c α s . Indeed, there have been several interesting studies of classical Yang-Mills fields under longitudinal expansion in the past years [2][3][4][5][6]. On the other hand, once the system has reached a local thermal equilibrium, or at least approximately isotropized [7], the matter in the mid-rapidity region (at sufficiently low p T ) is described by relativistic fluid dynamics [8,9]. There has been a sizeable and very successful program of numerical simulations of relativistic hydrodynamics.
It is well known that classical Yang-Mills theory can not reach thermalization due to the Rayleigh-Jeans catastrophe nor can it isotropize when exposed to rapid longitudinal expansion as noticed by [10]. Instead classical evolution drives the system further away from equilibrium, making the system less occupied (i.e. weaker fields) but more anisotropic, which has also been observed in the simulations of [4,5]. Therefore, there is a missing link between the physics of saturated gluon fields and fluid dynamics that needs to be bridged in order to fully carry the predictions of saturation physics to the hydrodynamic regime. This gap can be filled with the systematically improvable framework of effective kinetic theory (EKT) of [11].
The EKT faithfully describes to leading order in λf systems where the typical occupancies of gluons are not nonperturbative f ≪ 1/λ and have momentum significantly larger than the in-medium screening scale p 2 > m 2 ≡ λ p f (p)/p. At weak coupling, these conditions are certainly fulfilled in thermal equilibrium and there-fore we can describe the system all the way to the equilibrium. While these conditions are not fulfilled at the very earliest times, both EKT and classical Yang-Mills give an equally valid leading order description for a wide range of large but perturbative occupancies 1 ≪ f ≪ 1/λ [12,13]. Therefore a possible strategy to simulate the system through all time scales is to start the simulation with a classical Yang-Mills simulation (for Q s τ ∼ 1 and f ∼ 1/λ) and subsequently pass the system to EKT at some arbitrary time τ EKT Q s ≫ 1. Then, once EKT has brought the system sufficiently close to the thermal equilibrium, both hydrodynamics and EKT should give equivalent descriptions and the system can be passed to a hydrodynamical simulation at some arbitrary time τ hydro .
The proof of principle of such a procedure was shown in a series of papers [13][14][15] in an isotropic setting in the absence of expansion. Here, we present first results of numerical simulations of the EKT in a boost invariant 2+1D setting and demonstrate the connection with both classical Yang-Mills simulations and viscous hydrodynamics. We make a first attempt to model the early stages of heavy-ion collision starting from the overoccupied region all the way to the thermal equilibrium and find that after time Q s τ ∼ 5.0 the time evolution of energy density is described by viscous hydrodynamics to a better than 10% accuracy for a realistic coupling λ = 10 corresponding to α s ≈ 0.3.
The current result is, however, incomplete in two ways. Firstly, there have been significant steps to describe the melting of the strong color-fields in classical statistical field theory [2,3,6], but currently we do not have a reliable initial state from a 3+1D simulation to enter into our EKT simulation. Therefore we initialize our system at Q s τ ∼ 1 with the energy density given by 2+1D simulation [28] and vary the parameters of the initial condition to quantify the ignorance of the very early time dynamics. Secondly, certain non-perturbative chromo-Weibel instabilities [16] arising from anisotropic screening play a parametrically leading order role during the whole non-equilibrium evolution [17,18]. However, numerical simulations of classical fields [6] have not been able to see them beyond the very early times, suggesting that even if parametrically leading order, their effect is numerically small for values of λ that are phenomenologically interesting. Therefore, as a formally correct treatment of anisotropic screening is rather complicated [19][20][21][22] and not a fully solved problem, we in this first work will treat the screening in an isotropic way. Therefore, our results are correct to leading order in λ for systems that are close to isotropy (late times) and for large anisotropies our results have leading logarithmic accuracy apart from the instabilities.

METHODOLOGY
The EKT of [11] is defined though the effective Boltzmann equation for the color and spin averaged distribution function of gluons, and to leading order in λ it contains effective 2 ↔ 2 scattering and 1 ↔ 2 splitting terms. (1) We restrict ourselves to azimuthally symmetric distributions but allow anisotropy in the z-direction so that it is enough to specify [23]. The 2 ↔ 2 effective scattering term reads where ν = 2d A and g p ≡ 1 + f p . dΓ P S is the integral measure over the phase space of 2 ↔ 2 processes, singling out the z-direction In these coordinates the angles of incoming and outgoing momenta read This relation holds for x k , x p ′ , and x k ′ with cos φ p ′ q = cos φ pq , cos φ k ′ q = cos φ kq , and t ≡ ω 2 − q 2 . The effective matrix element for most kinematics is the ordinary vacuum tree-level element For small t ∼ m 2 or u ∼ m 2 the tree-level element has an IR divergence that is regulated by the physics of screening. In the formulation of [11], the screening is treated by replacing the naive tree-level matrix element by the retarded Hard Thermal Loop (HTL) resummed expression in the soft kinematic region. The resummed propagator in the anisotropic case, however, contains a non-integrable singularity, signalling the presence of an instability in the system. Therefore, in the anisotropic system this procedure is not fully satisfactory. Here, in this work we will restrict ourselves to isotropic screening and use the prescription derived in [13], replacing in the denominator (similarly for u) with ξ 0 = e 5/6 / √ 8, which is leading order accurate in the case of isotropic distribution.
The effective 1 ↔ 2 collinear splitting term reads where p ′ = p − k ′ , and the effective splitting rate reads Again, consistently assuming only isotropic screening, the function F is given by the solution to the linear integral equation We solve this equation using an efficient numerical method introduced in [24].
For the numerical solution of the EKT, we discretize n xp,p = 4πp 2 /(2π) 3 f xp,p on a 2D grid and Monte Carlo estimate the integrals of Eqs. (3,9). We use a logarithmically spaced grid for both p and x p with at least 250 grid points in angular direction and at least 100 in the p-direction. We have varied the number of grid points to verify that the results are insensitive to the number of grid points. Our algorithm, the discrete-p method of [13], conserves energy exactly and has exactly correct particle number violation for 1 ↔ 2 term.  Fig. 2. The simulations at finite coupling reach thermal equilibrium located at points indicated by the black crosses.

RESULTS
We shall now apply EKT to simulate the prethermal evolution of the expanding fireball created in a heavyion collision. In saturation framework, the initial condition is typically described in terms of "gluon liberation coefficient" c and mean transverse momentum p T /Q s [26,27]. The gluon liberation coefficient is proportional to the total gluon multiplicity per unit rapidity after the classical fields have decohered and can be described in terms of quasi-particles. Lappi [28] finds in JIMWLK evolved MV model values relevant for heavyion collisions relevant for LHC of roughly p T ≈ 1.8Q s and c ≈ 1.25 extracted at time Q s τ = 12 from a 2D classical Yang-Mills simulation. By construction the distribution then has p z = 0. But it is has been noted [25] that certain plasma instabilities will broaden the distribution in p z in a time scale Qτ ∼ 1/ log 2 (λ −1 ). Therefore, as a rough estimate of the initial condition we instead take somewhat arbitrarily our initial condition at the time Qτ = 1 to be choosing A such that comoving energy density τ ǫ = p T dN/d 2 xdy is fixed. We then vary ξ = 4, 10 to quantify our ignorance of the initial nonperturbative dynamics. Figure 1 displays a set of trajectories from simulations with varying λ and ξ = 4, 10 on a plane of mean occupancy (weighted by the energy of particles) and anisotropy measured by the ratio of the transverse and longitudinal pressures P T /P L . The line with λ = 0 corresponds to the classical field limit λ → 0 with fixed λf , which is obtained in EKT by including only the highest power of f 's in Eqs. (2,9). The classical field theory can not thermalize and indeed instead it flows to a stationary scaling solution. By performing classical Yang-Mills simulations Berges et. al have established that the scaling solution can described by a scaling form of the distribution function [4], where f S is approximately constant as a function of time. This behavior is demonstrated in Figure 2 where we plot a section of rescaled distribution function f S measured at various times as a function forp z ≡ (Q s τ ) 1/3 p z at fixed p ⊥ following Berges et al.. Our results corroborate that such a scaling solution exists at late times within the classical approximation and we observe that the scaling regime is reached after a time Q s τ ∼ 15. Moving on to the finite but small coupling λ = 1, 0.5, we see qualitative agreement with the parametric picture of bottom-up thermalization of [10]: there are three distinct stages of evolution visible. In the first stage the classical evolution drives the system more anisotrpic and less occupied. Once the occupancies reach f ∼ 1, there is a qualitative change in the dynamics of the system as Bose enhancement is lost. This has the effect that ansiotropy freezes but the system still continues to get more dilute. Only in the last stage which is characterized by a radiational break up of the particles at the scale Q s , does the trajectory turn back and reach thermal equilibrium, denoted by the black crosses in the Figure 1.
For larger values of coupling λ = 5.0, 10, these features become however less pronounced and the system takes more straight trajectory towards equilibrium. It should be noted that for these values of λ the assumption of p ≫ m is not satisfied and large NLO corrections are to be expected.

Approach to hydrodynamics
We expect that late time evolution should be described by relativistic hydrodynamics. Under flow with translational invariance along transverse directions and boost invariance, the hydrodynamical relations read to second order in gradients [29] ∂ τ ǫ = − 4 3 with longitudinal and transverse pressure P L = 1 3 ǫ − Φ and P T = 1 3 ǫ + 1 2 Φ. First order hydrodynamics corresponds to setting Φ = 4η 3τ in Eq. (15). At weak coupling, the transport coefficients η, τ Π and λ 1 are known [30,31] leaving zero free parameters to fit. besides a time when the hydrodynamics is initialized. We initialize the 1st order hydrodynamics at the latest time we have in our simulation and integrate Eq. (15) backwards in time. For 2nd order hydrodynamics integrating backwards is highly unstable and we initialize the energy density at some arbitrary earlier time and integrate forwards in time.
In Figure 3 we examine the validity of the hydrodynamical expansion at small λ = 1 and at realistic intermediate λ = 10 (α s ≈ 0.3) values of coupling. In both cases we see that the evolution of the components of the energy momentum tensor asymptotes to their hydrodynamical values. In case of λ = 1, the hydrodynamical behaviour is reached only at a rather late time Q s τ ∼ 2000. We have checked that including 2nd order terms before this time does not make the convergence significantly better; before this time the evolution differs qualitatively from the hydrodynamical prediction. However, rather remarkably, for λ = 10 even 1st order hydrodynamics gives a very good description of the data all the way to very early times Q s τ ∼ 10 (corresponding to τ ∼1fm/c for Q s = 2.0GeV) where the ratio of the pressures is still as large as P T /P L ≈ 5. In addition, including the second order terms significantly improves the convergence. Indeed, we find that initializing the 2nd order hydrodynamics at Q s τ = 1 leads to only 10% relative error in the energy density at late times.

DISCUSSIONS AND CONCLUSION
The parametric estimate of Baier et al. [10] for the time when the hydrodynamic behaviour sets is Q s τ ∼   λ −13/5 . This arises from equating the age of the system τ with the time scale τ Q it takes to affect appreciably the scale Q s in a thermal bath whose temperature depends on this time T (τ ) ∼ λ −1/4 Q s (Q s τ ) −1/4 according to conservation of comoving energy density. In [10] it was assumed that the rate for affecting the scale Q s is LPM suppressed [32], giving τ Q ∼ (Q s /T ) 1/2 /λ 2 T . A selfconsistent solution of these equations gives the aforementioned estimate of [10]. Arnold and Lenaghan [33] noted that as the scattering rate in thermal plasma is τ Q ∼ 1/λ 2 T , there can be no process that would make the estimate faster than Q s τ ∼ λ −7/3 . We have examined the validity of these scaling estimates by plotting the difference of the energy density obtained from the simulation and the 1st order hydrodynamic estimate and find that both of these estimates describe the data poorly. However, if we assume that the approach is governed by the hydrodynamical relaxation time τ Q ∼ τ Π ∼ η/sT we get an estimate τ ∼ λ 1/3 (η/s) 4/3 /Q s . Figure 4 displays the deviation of the hydrodynamics as a function of rescaled time. In particular for intermediate couplings λ = 5, 10 the overlap of the different curves indicates correct scaling. Note that this estimate is parametrically the same as the estimate of Arnold and Lenaghan as parametrically η/s ∼ λ −2 . However, there are large corrections beyond the parametric estimate in η/s due to which it is important to use the full numerical result instead of the simpler parametric estimate. We however believe (in the absence of plasma instabilities) that the correct scaling at sufficiently small λ is that of [10]. This estimate is however based on a large scale separation of T (τ hydro )/Q s ∼ λ 2/5 . Numerically this ratio is not very large even in our simulation with small λ = 1.0, where the ratio is T /Q s ≈ 0.15 (at Q s τ = 1800) as can be inferred from Figure 3. Therefore we believe that the scaling predicted by [10] sets in only at significantly smaller couplings.
In conclusion, we have shown how a far from equilibrium overoccupied configuration of gluons reaches hydrodynamical flow under longitudinal expansion in a weak coupling setting that is systematically improvable. It has been demonstrated by Chesler and Yaffe within AdS/CFT framework that at large values of 't Hooft coupling, hydrodynamics is surprisingly robust even in the presence of large anisotropy [34]. The main deliverable of this paper is to show that this robustness is present also in the weak coupling picture extrapolated to intermediate couplings relevant for heavy-ion collisions.