Elliptic Flow from Non-equilibrium Initial Condition with a Saturation Scale

A current goal of relativistic heavy ion collisions experiments is the search for a Color Glass Condensate as the limiting state of QCD matter at very high density. In viscous hydrodynamics simulations, a standard Glauber initial condition leads to estimate $4\pi \eta/s \sim 1$, while a Color Glass Condensate modeling leads to at least a factor of 2 larger $\eta/s$. Within a kinetic theory approach based on a relativistic Boltzmann-like transport simulation, we point out that the out-of-equilibrium initial distribution proper of a Color Glass Condensate reduces the efficiency in building-up the elliptic flow. Our main result at RHIC energy is that the available data on $v_2$ are in agreement with a $4\pi \eta/s \sim 1$ also for Color Glass Condensate initial conditions, opening the possibility to describe self-consistently also higher order flow, otherwise significantly underestimated, and to pursue further the search for signatures of the Color Glass Condensate.

Ultra-relativistic heavy-ion collisions (uRHICs) at the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) create a hot and dense system of strongly interacting matter. In the last decade it has been reached a general consensus that such a state of matter is not of hadronic nature and there are several signatures that it is a strongly interacting quark-gluon plasma (QGP) [1][2][3]. A main discovery has been that the QGP has a very small shear viscosity to density entropy, η/s, which is more than one order of magnitude smaller than the one of water [4,5], and close to the lower bound of 1/4π conjectured for systems at infinite strong coupling [6]. A key observable to reach such a conclusion is the so-called elliptic flow, v 2 = cos(2ϕ p ) = (p 2 x − p 2 y )/(p 2 x + p 2 y ) , with ϕ p being the azimuthal angle in the transverse plane and the average meant over the particle distribution. In fact, the expansion of the created matter generates a large anisotropy of the emitted particles that can be primarily measured by v 2 . Its origin is the initial spatial eccentricity, ǫ x = x 2 − y 2 / x 2 + y 2 , of the overlap region in non-central collisions. The observed large v 2 is considered a signal of a very small η/s because it means that the system is very efficient in converting ǫ x into an anisotropy in the momentum space v 2 , a mechanism that would be strongly damped in a system highly viscous that dissipates and smooths anisotropies [7,9,11]. Quantitatively both viscous hydrodynamics [7-9, 12, 13, 15], and transport Boltzmannlike approaches [16][17][18][19][20] agree in indicating an average η/s of the QGP lying in the range 4πη/s ∼ 1 − 3.
The uRHIC program offers the tantalizing opportunity to explore the existence of an exotic state, namely the Color Glass Condensate (CGC) [21,22], see [23,24] for reviews. Such a state of matter would be primarily generated by the very high density of the gluon parton distribution function at low x (parton momentum fraction), which triggers a saturation of the gluon distribution func-tion at a p T below the saturation scale, Q s [25]. Even if at first sight surprisingly, the study of the shear viscosity η/s of the QGP and the search for the CGC are related. In fact, the main source of uncertainty for η/s comes from the unknown initial conditions of the created matter [8] and confirmed later by further works [10,12,26].
A simple geometrical description through the Glauber model [28] predicts a ǫ x smaller at least 25-30% than the eccentricity of the CGC, for most of the centralities of the collisions, see for example results within the Kharzeev-Levin-Nardi (KLN) model [27,29,30], factorized KLN (fKLN) model [31], Monte Carlo KLN (MC-KLN) model [31,32] and dipole model [27,33]. The uncertainty in the initial condition translates into an uncertainty on η/s of at a least a factor of two as estimated by mean of several viscous hydrodynamical approaches [8,10,12,26]. More explicitely, the experimental data of v 2 (p T ) at the highest RHIC energy are in agreement with a fluid at 4πη/s ∼ 1 according to viscous hydrodynamics simulation, assuming a standard Glauber initial condition. Assuming an initial fKLN or MC-KLN space distribution the comparison favors a fluid at 4πη/s ∼ 2. The reason is the larger initial ǫ x of the fKLN, which leads to larger v 2 unless a large η/s is considered. However, in [10,26] it has been shown that viscous hydrodynamics fails to reproduce both v 2 and v 3 = cos(3ϕ p ) if the same 4πη/s ∼ 2 is assumed. At variance a Glauber initial condition can account for both with the same 4πη/s = 1. However the indirect effect on ǫ x is not a unique and solid prediction of the CGC modellings; for example, the approach based on the solution of the Classical Yang-Mills (CYM) equations predict a somewhat smaller initial eccentricity [34,36]. Very recently employing a x − space distribution inspired to the CYM approach in a viscous hydrodynamical approach [35] it has been show that not only v 2 but also higher harmonics can be correctly predicted with a 4πη/s ∼ 1.5 instead of ∼ 2, which is in qualitative agreement with the fact that CYM tend to predict quite smaller ǫ x with respect to fKLN. However our present studies focus on the effect of the initial nonequilibrium in p−space, an issue discarded in all previous studies including the recent ones [35,36].
In this Letter, we point out that the implementation of the melted CGC in hydrodynamics takes into account only the different space distribution with respect to a geometric Glauber model, discarding the key and more peculiar feature of the damping of the distribution for p T below the Q s saturation scale. We have found by mean of kinetic theory that this has a pivotal role on the build-up of v 2 . We adopt the model which was firstly introduced by Kharzeev, Levin and Nardi [29] (KLN model) even if in the regime of over saturation in A + A collisions some aspects are better caught by a CYM approach [48]. In particular, to prepare the initial conditions of our simulations we refer to the factorized-KLN (fKLN) approach as introduced in [31,32]. This will allow for a direct comparison with viscous hydrodynamics results. in which the coordinate space distribution function of gluons arising from the melted CGC is assumed to be where Φ corresponds to the momentum space distribution in the k T factorization hypothesis [37,38], Here x 1,2 = p T exp(±y)/ √ s and the ultraviolet cutoff p T = 3 GeV/c assumed in the p T integral in Eq. (1); α S denotes the strong coupling constant, which is computed at the scale Q 2 = max(k 2 T , (p T − k T ) 2 ) according to the one-loop β function but frozen at α s = 0.5 in the infrared region as in [30,33,39]. In Eq. (1) p A,B denote the probability to find one nucleon at a given transverse coordinate, where σ in is the inelastic cross section and T A corresponds to the usual thickness function of the Glauber model.
The main ingredient to specify in Eq. (2) is the unintegrated gluon distribution function (uGDF) for partons coming from nucleus A, which is assumed to be: where we see the peculiar feature of the CGC that is the saturation of the distribution for p T < Q s ; a similar equation holds for partons belonging to nucleus B. Following [31] we take the saturation scale for the nucleus A as with λ = 0.28, and similarly nucleus B. This choice is the one adopted in fKLN or MC-KLN and in hydro simulations [12,31] to study the dependence of v 2 (p T ) on η/s. Using Eqs. (4) and (1) we find that Q s ≈ 1.4 GeV where the average is understood in the transverse plane. We employ transport theory as a base of a simulation code of the fireball expansion created in relativistic heavy-ion collision [19,20,40,41], therefore the time evolution of the gluons distribution function f (x, p, t) evolves according to the Boltzmann equation: where d 3 p k = 2E k (2π) 3 dΓ k and M corresponds to the transition amplitude.
At variance with the standard use of transport theory, we have developed an approach that, instead of focusing on specific microscopic calculations or modelings for the scattering matrix, fixes the total cross section in order to have the wanted η/s. In Ref. [42] it has been shown in 1+1D such an approach is able to recover the Israel-Stewart viscous hydrodynamical evolution when η/s is sufficiently small. In 3+1D some of the authors has studied the analytical relation between η, temperature, cross section and density and as shown in [41,44], the Chapmann-Enskog approximation supplies such a relation with quite good approximation [43], in agreement with the results obtained using the Green Kubo formula. Therefore, we fix η/s and compute the pertinent total cross section by mean of the relation which is valid for a generic differential cross section dσ/dt ∼ α 2 s /(t − m 2 D ) 2 as proved in [44]. In the above equation a = T /m D , with m D the screening mass regulating the angular dependence of the cross section, while g(a) = 1 50 with K n the Bessel function and h corresponding to the ratio of the transport and the total cross section. The maximum value of g, namely g(m D → ∞) = h(m D → ∞) = 2/3, is reached for isotropic cross section; a smaller value of g(a) means that a higher σ tot is needed to reproduce the same value of η/s. However, we notice that in the regime were viscous hydrodynamic applies (not too large η/s and p T ) the specific microscopic detail of the cross section is irrelevant and our approach is an effective way to employ transport theory to simulate a fluid at a given η/s. From the operative point of view, keeping η/s constant in our simulations is achieved by evaluating locally in space and time the strength of the cross section by means of Eq. (6), where both parton densities and temperature are computed locally in each cell. To realize a realistic freeze-out, when the local energy density reaches the cross-over region, the η/s increases linearly to match the estimated hadronic viscosity, as described in [20,41]; this affects in the same way all the cases considered in the following. In the following, we will consider three different types of initial distribution function in the phase-space, two of which are the one employed till now for the investigation of the η/s, while the third one is the genuine novelty of the present study. Furthermore, we refer to Au + Au collision at √ s = 200 AGeV and b = 7.5 fm. In this case, our result for initial eccentricity in the fKLN model is ǫ x = 0.357 (which is in agreement with MC-KLN [31] result used in hydro simulations). The standard initial condition for simulations of the plasma fireball created at RHIC is a x-space distribution given by the Glauber model and a p-space thermalized spectrum in the transverse plane at a time τ = 0.6 fm/c with a maximum initial temperature T 0 = 340MeV. In this case, for a standard mixture of N part and N coll we find ǫ x = 0.284. We will refer to this case as Th-Glauber. Instead the study of the impact of an initial CGC state has been performed considering an x-space distribution given by the fKLN (or MC-KLN), while in the momentum space the spectrum has been considered thermalized at τ 0 ∼ 0.6fm/c; we refer to this case as Th-fKLN and it is represents the case implemented in hydrodynamics, that has lead to the conclusion that the CGC suggests a 4πη/s ∼ 2 [8,10,12,26]. The third initial conditions is the full fKLN initial conditions where, beyond the x-space, the saturated distribution in p-space is implemented as well, see Fig. 1 solid thick line. As initial time we take τ 0 = 0.15fm/c because in this case there is no pre-assumption of thermalization. This is not usually considered in hydrodynamics because there it is implicitly assumed a distribution function in p-space in local equilibrium, at least in the transverse plane. The choice of τ 0 for the glasma-like initial condition is inspired by the recent results of [50][51][52] where it is discussed that, even if at τ = 0 + the longitudinal pressure of the glasma is negative, within a time τ ≈ 1/Q s it becomes positive thus making a description of the expanding system in terms of a partonic distribution function quite reliable. For all the previous cases, as usually done, a Bjorken scaling at the initial time is assumed, identifying momentum rapidity y z with space rapidity η s = arctg −1 (z/τ ). For all the case considered the multiplicity dN/dη at mid rapidity has been fixed initially equal to 400 to approximately match the experimental data that for the impact parameter considered corresponds to about N part = 150 in [45].
In Fig. 1 we plot the initial spectra for the fKLN (thick solid line), Th-fKLN (dashed line) and Th-Glauber (thin solid line) at their respective initial times τ 0 , and the spectrum of the fKLN model after a time evolution ∆τ = 0.8 fm/c (dashed green line). We notice that initially the fKLN spectrum is quite far from a thermalized spectrum; in fact, it embeds the saturation effects which are proper of the melted CGC. Neverthelsess, the spectrum evolves to a thermalized one within 1 fm/c. Such a feature is confirmed by the inset of Fig. 1, where the quantity T * · τ 1/3 is shown, with T * = E/3N representing, the temperature in the case of a thermalized system. It is known that in the case of 1D expansion a thermalized system should keep such a ratio constant. We find that in the case of the fKLN (solid green line) the product T * ·τ 1/3 is strongly dependent on time because the system is quite far from equilibrium; however at τ ∼ 0.8 fm/c both the value and the time evolution become indistinguishable from the thermal cases represented by the Th-Glauber and Th-fKLN. We notice a little adjustment also for these cases that we have indicated as thermal. The reason is that the initial spectra are thermal only in the transverse plane, but there is a boost invariance in the longitudinal direction. This causes a little re-adjustment that would disappear assuming a thermal spectral also in the longitudinal direction.
Our results on thermalization time are not in disagreement with earlier studies showing that two-body collisions are insufficient to achieve a fast thermal equilibrium [46]. In fact in that case a perturbative QCD two-body cross section is employed, corresponding to η/s about one order of magnitude larger than in our case [44], while here we normalize the cross section to get an η/s = 0.08 that implies scattering rates very large that induce a fast thermalization.
We have studied also the isotropization of the pressures of the expanding system. In order to do this we have computed the relevant components of the energymomentum tensor, where f corresponds to the invariant distribution function. In our simulations the energy-momentum tensor is defined in each cell; we then define transverse and longitudinal pressures, P T and P L respectively, as where the integration is restricted to the region Ω defined by |x| < 2.45 fm, |y| < 4.55 fm and |η| < 0.5, and V is the volume of such a region. In the lower panel of Fig. 2 we plot our results about P L /P T as a function of time, for the three initial conditions considered till now. As we have already explained, the initial time for the simulation depends on the particular initial condition used. Our findings suggest that independently on the initial condition implemented, the system becomes isotropic within 1 fm/c. A similar result has been obtained recently using both 2 ↔ 2 and 2 ↔ 3 collisions with pQCD cross sections and α s = 0.6 which should correspond to η/s ≈ 0.1 [49]. This also mean that the main effect we observe on the elliptic flow is not driven by the isotropization itself. This result is not in disagreement with other studies [50,51] which show how the expanding glasma becomes almost isotropic within few fm/c if the coupling is strong enough, or in general the dynamics sufficiently nondissipative [52].
In the left panel of Fig. 3 we plot the v 2 (p T ) for the case of Th-Glauber (thin solid line) and Th-fKLN (thick solid line) at a fixed 4πη/s = 1. The Glauber initial condition reproduces quite well the data (circles); in the case of Th-fKLN (thick solid line) one gets a too large v 2 and for such initial conditions the agreement with the data is achieved only if the η/s is increased by a factor of two (dashed line). These results are in agreement with the ones obtained from viscous hydrodynamics [8,10,12,26], showing the solidity and consistency of our transport approach at fixed η/s.
In the right panel of Fig. 3 we present our novel result for the fKLN model, when the CGC distribution function is implemented in both the x and p spaces. We find that fKLN with a 4πη/s = 1 (thick solid line) gives a v 2 (p T ) quite similar to the Th-Glauber, while in such a case if 4πη/s = 2, dashed line, the v 2 (p T ) would be too small. Our interpretation is that the initial larger ǫ x is compensated by the key feature of an almost saturated initial distribution in p-space below the saturation scale Q s . In other words the initial out-of-equilibrium fKLN distribution reduces the efficiency in converting ǫ x into v 2 . In fact, the elliptic flow can be understood as a larger slope of the momentum spectrum in the out of plane x direction with respect to the y one caused by a larger pressure in the x direction due the elliptical shape. The net effect in terms of the difference of the particle yields between the two directions is larger if the spectra are decreasing exponentially with respect to the case in which they are nearly flat as a function of p T . A detailed study is in preparation, but the result we present in this Letter shows that the initial out-of-equilibrium function implied by the CGC cannot be discarded in studying the build of the v 2 and the effect change significantly the estimate of η/s. We however notice that in the context of viscous hydrodynamics it is possible to include initial non-equilibrium conditions by introducing an initial non-vanishing value for the viscous tensor Π µν (t 0 ). This has been seen to have a quite small impact on the v 2 at least on the bulk of the system [8,14,15]. However to our knowledge it has never been investigated the relation, if any, between the non-equilibrium implied by Π µν (t 0 ) and the one implied by the saturation scale in the distribution function f (x, p) in CGC inspired models.
In Fig. 4, we show the build-up of v 2 at a fixed p T = 1.5 GeV as a function of time, τ − τ 0 , for all the initial conditions considered. We find that the rate of increase of v 2 for the thermal distributions is large from the very beginning of the evolution, while it is quite reduced for fKLN. After about 1−1.5fm/c, roughly corresponding to the thermalization also for the full fKLN distribution, the increase of v 2 with time also for fKLN becomes very similar to Th-fKLN. We notice that when plotted as a function of τ − τ 0 the evolution of the space eccentricity is quite similar between T h − f KLN and f KLN , hence the differences observed are mainly driven by the p-space distributions. This observation further confirms that it is the initial out-of -equilibrium and nearly saturated distribution that dampens the efficiency in converting the space eccentricity into the v 2 . Therefore, even if the thermalization sets in quite quickly, as assumed in hydrodynamics τ ∼ 0.8 fm/c, such a time duration cannot be neglected in studying the impact of initial conditions.
In conclusion, our study based on kinetic theory shows that the elliptic flow in a system depends not only on the pressure gradients and the η/s of the system, but also on the initial distribution in momentum space. A p-distribution with a saturation behavior generates smaller v 2 with respect to the thermal one. This result is quite general, and we expect it should be valid, besides QGP in uRHICs, for systems like cold atoms in a magnetic trap which are characterized by a value of η/s close to the QGP one [47]. In the specific case of the KLN matter studied here, the effect of the initial non-equilibrium distribution affects the estimate of η/s of about a factor of two. The relevance of our results is further enhanced by the fact that Th-fKLN with 4πη/s ∼ 2 would generate a low v 3 with respect to the available data, which is the main conclusion of [26]. We notice that in the present kinetic approach the quantum nature of gluons has been discarded. This is justified at RHIC because in this Letter we have focused on an effect which is dominant at p T > 0.5 GeV, where the f (x, p) is still smaller than unity. At LHC, or anyway at small p T , it would be necessary to include (1 + f ) terms in Eq.(5) that could drive the system toward a Bose-Einstein condensate [53].