Fluctuating Glasma initial conditions and flow in heavy ion collisions

Fluctuating Glasma initial conditions and flow in heavy ion collisions


We compute initial conditions in heavy-ion collisions within the Color Glass Condensate (CGC) framework by combining the impact parameter dependent saturation model (IP-Sat) with the classical Yang-Mills description of initial Glasma fields. In addition to fluctuations of nucleon positions, this IP-Glasma description includes quantum fluctuations of color charges on the length-scale determined by the inverse nuclear saturation scale . The model naturally produces initial energy fluctuations that are described by a negative binomial distribution. The ratio of triangularity to eccentricity is close to that in a model tuned to reproduce experimental flow data. We compare transverse momentum spectra and of pions from different models of initial conditions using relativistic viscous hydrodynamic evolution.

A large uncertainty in the hydrodynamical description of ultrarelativistic heavy ion collisions is our imperfect knowledge of multi-gluon states in the nuclear wavefunctions and the early-time dynamics of gluon fields after the collision. In heavy ion collisions, studies of observables sensitive to harmonics of hydrodynamic flow distributions provide constraints both on the low shear viscosity of the Quark-Gluon Plasma (QGP) and the initial state dynamics Alver and Roland (2010); Alver et al. (2010); Petersen et al. (2010); Schenke et al. (2011a, b); Qiu and Heinz (2011); ?; Bleicher et al. (2011); Schenke et al. (2011c). This situation is analogous to studies of the cosmic microwave background Komatsu et al. (2009), wherein inhomogeneities in the observed power spectrum are sensitive to primordial quantum fluctuations.

An ab initio framework for multi-particle production is the Color Glass Condensate (CGC) Gelis et al. (2010) wherein the initial state dynamics is described by flowing Glasma gluon fields Kovner et al. (1995a); ?; Krasnitz and Venugopalan (2000); ?. There are several sources of quantum fluctuations that can influence hydrodynamic flow on an event-by-event basis. An important source of fluctuations, generic to all models of quantum fluctuations, are fluctuations in the distributions of nucleons in the nuclear wavefunctions. In addition there are fluctuations in the color charge distributions inside a nucleon. This, combined with Lorentz contraction, results in “lumpy” transverse projections of color charge configurations that vary event to event. The scale of this lumpiness is given on average by the nuclear saturation scale which corresponds to distance scales smaller than the nucleon size Kowalski et al. (2008). For each such configuration of color charges, the Quantum Chromo-Dynamics (QCD) parton model predicts dynamical event-by-event fluctuations in the multiplicities, the impact parameters and the rapidities of produced gluons Miettinen and Pumplin (1978).

All these sources of fluctuations are captured in the CGC Glasma flux tube picture. The relevant feature of this scenario is that long range rapidity correlations from the initial state wavefunctions, coherent over in the transverse plane, are efficiently transmitted into hydrodynamic flow of the final state quark-gluon matter Voloshin (2006); Shuryak (2007).

Recently, Monte-Carlo Glauber-type models (MC-Glauber) and Monte-Carlo implementations of the Kharzeev-Levin-Nardi-model (MC-KLN) Hirano et al. (2006); Drescher et al. (2006); ? have been compared to experimental data on elliptic and triangular moments of the flow distribution. While both types of models treat fluctuations in nucleon positions identically, the Glauber model implementations do not specify a mechanism for multi-particle production which would constrain the initial energy density distribution. MC-Glauber initial conditions Esumi (2011); Qiu and Heinz (2011); ? can be tuned to reproduce data on both elliptic and triangular flow from RHIC and the LHC. The MC-KLN model is motivated by the CGC with approximations that will be discussed further below. It requires larger values of the viscosity to entropy density ratio () relative to the Glauber model values to describe elliptic flow data. This however leads it to underpredict triangular flow data.

Odd flow harmonics are entirely driven by fluctuations; it is therefore essential to have a realistic description of quantum fluctuations in multi-particle production to properly describe the final state dynamics. Towards this end, we will consider in this letter the impact parameter dependent saturation model (IP-Sat) Bartels et al. (2002); Kowalski and Teaney (2003) to determine fluctuating configurations of color charges in two incoming highly energetic nuclei. This model is formally similar to the classical CGC McLerran-Venugopalan (MV) model of nuclear wavefunctions McLerran and Venugopalan (1994a); ?; ?, but additionally includes Bjorken and impact parameter dependence through eikonalized gluon distributions of the proton that are constrained1 by HERA inclusive and diffractive e+p deeply inelastic scattering (DIS) data Kowalski et al. (2006). Most importantly, the model is in excellent agreement with data on n-particle multiplicity distributions in p+p collisions at RHIC and the LHC and in A+A collisions at RHIC Tribedy and Venugopalan (2011); ?, an essential requirement for microscopic models. The MC-KLN model does not contain these features; a scheme to introduce fluctuations in the model has only been discussed recently Dumitru and Nara (2012).

The color charges in the IP-Sat/MV model act as local sources for small classical gluon Glasma fields. These are determined by solving the classical Yang-Mills (CYM) equations , with the color current generated by a nucleus A (B) moving along the () direction.2 The solution in light cone gauge are the pure gauge fields  McLerran and Venugopalan (1994a); ?; ?; Jalilian-Marian et al. (1997); Kovchegov (1996) and . Here is a path ordered Wilson line in the fundamental representation, where the infrared regulator (of order ) incorporates color confinement at the nucleon level.

The initial condition for a heavy-ion collision at time is given by the solution of the CYM equations in Schwinger gauge , a natural choice because it interpolates between the light cone gauge conditions of the incoming nuclei. It has a simple expression in terms of the gauge fields of the colliding nuclei Kovner et al. (1995b, a):


and3 , . In the limit , , with the longitudinal component of the electric field. At , the only non-zero components of the field strength tensor are the longitudinal magnetic and electric fields, which can be computed non-perturbatively. They determine the energy density of the Glasma at at each transverse position in a single event Krasnitz and Venugopalan (2000); ?; Lappi (2003). The Glasma distribution computed from the CYM equations4 (IP-Glasma) is matched event-by-event to viscous relativistic hydrodynamics Schenke et al. (2011a, b) to compute harmonics of the flow distributions.

We will now discuss details of the computation. Nucleon positions in the nucleus are sampled from a Fermi distribution. The saturation scale is determined from the IP-Sat dipole cross section for each nucleon, where is the impact parameter relative to each nucleon’s center. The color charge squared per unit transverse area is proportional5 to . For the nucleus, is obtained Krasnitz et al. (2003) by adding the individual nucleon at the same and transverse position in the nucleus.

Figure 1: The IP-Glasma event-by-event distribution in energy for fm on the lattice compared to different functional forms. The negative binomial distribution (NBD) gives the best fit.

The lattice formulation of the Glasma initial conditions in Eq. (1) was first given in Krasnitz and Venugopalan (1999). On a transverse lattice, random color charges6 are sampled from


where the indices7 represent a discretized coordinate Lappi (2008). For the large nuclei we consider the use of such local Gaussian color charge distributions is a valid approximation8. The path ordered Wilson line is discretized as


To each lattice site we assign two SU() matrices and , each of which defines a pure gauge configuration with the link variables where indicates a shift from by one lattice site in the transverse direction. The link variables in the future lightcone , are determined from solutions of the lattice CYM equations at ,


where are the generators of in the fundamental representation (The cell index is omitted here). The equations (4) are highly non-linear and for are solved iteratively.

The total energy density on the lattice at is given by


where the first term is the longitudinal magnetic energy, with the plaquette given by . The explicit lattice expression for the longitudinal electric field in the second term can be found in Refs. Krasnitz and Venugopalan (1999); Romatschke and Venugopalan (2006). We note that the boost-invariant CYM framework neglects fluctuations in the rapidity direction. Anisotropic flow at mid-rapdity is dominated by fluctuations in the transverse plane but fluctuations in rapidity could have an effect on the dissipative evolution; the framework to describe these effects has been developed Dusling et al. (2011) and will be addressed in future work. Other rapidity dependent initial conditions are discussed in Ref. Magas et al. (2001); ?.

In Fig. 1 we show the event-by-event fluctuation in the initial energy per unit rapidity. The mean was adjusted to reproduce particle multiplicities after hydrodynamic evolution. This and all following results are for Au+Au collisions at RHIC energies () at midrapidity. The best fit is given by a negative binomial (NBD) distribution, as predicted in the Glasma flux tube framework Gelis et al. (2009); our result adds further confirmation to a previous non-perturbative study Lappi et al. (2010). The fact that the Glasma NBD distribution fits p+p multiplicity distributions over RHIC and LHC energies Tribedy and Venugopalan (2011); ? lends confidence that our picture includes fluctuations properly.

We now show the energy density distribution in the transverse plane in Fig. 2. We compare to the MC-KLN model and to an MC-Glauber model that was tuned to reproduce experimental data Schenke et al. (2011a, c). In the latter, for every participant nucleon, a Gaussian distributed energy density is added. Its parameters are the same for every nucleon in every event, with the width chosen to be to best describe anisotropic flow data. We will also present results for a model where the same Gaussians are assigned to each binary collision. The resulting initial energy densities differ significantly. In particular, fluctuations in the IP-Glasma occur on the length-scale , leading to finer structures in the initial energy density relative to the other models. As noted in Dumitru and Nara (2012), this feature of CGC physics is missing in the MC-KLN model.

Figure 2: (Color online) Initial energy density (arbitrary units) in the transverse plane in three different heavy-ion collision events: from top to bottom, IP-Glasma, MC-KLN and MC-Glauber Schenke et al. (2011c) models.

We next determine the participant ellipticity and triangularity of all models. Final flow of hadrons is to good approximation proportional to the respective Qin et al. (2010), which makes these eccentricities a good indicator of what to expect for . We define


where is the energy density weighted average. The results from averages over events for each point plotted are shown in Fig. 3. The ellipticity is largest in the MC-KLN model and smallest in the MC-Glauber model with participant scaling of the energy density (). The result of the present calculation lies in between, agreeing well with the MC-Glauber model using binary collision scaling (). We note however that this agreement is accidental; binary collision scaling of eccentricities, as shown explicitly in a previous work applying average CYM initial conditions Lappi and Venugopalan (2006), does not imply binary collision scaling of multiplicities.

The triangularities are very similar, with the MC-KLN result being below the other models for most impact parameters. Again, the present calculation is closest to the MC-Glauber model with binary collision scaling. There is no parameter dependence of eccentricities and triangularities in the IP-Glasma results shown in Fig. 3. It is reassuring that both are close to those from the MC-Glauber model because the latter is tuned to reproduce data even though it does not have dynamical QCD fluctuations.

We have checked that our results for are insensitive to the choice of the lattice spacing , despite a logarithmic ultraviolet divergence of the energy density at Lappi (2006). They are furthermore insensitive to the choice of , the ratio , and the uncertainty in Bjorken at a given energy.

Figure 3: (Color online) Average participant ellipticity (upper panel) and triangularity (lower panel) of the initial state. This calculation (circles), MC-KLN (squares), Glauber implementation with participant and binary collision scaling (triangles).

Finally, in Fig. 4 we present results for the transverse momentum spectrum and anisotropic flow of thermal pions after evolution using music Schenke et al. (2010, 2011a) with boost-invariant initial conditions and shear viscosity to entropy density ratio . Average maximal energy densities of all models were normalized to assure similar final multiplicities. More pronounced hot spots, as emphasized previously Gyulassy et al. (1997), affect the particle spectra obtained from flow, leading to harder momentum spectra in the present calculation compared to MC-KLN and MC-Glauber models. Differences in and are as expected from the initial eccentricities of the different models.

Figure 4: (Color online) Thermal transverse momentum spectra (upper) and anisotropic flow coefficients , , and as functions of (lower) from IP-Glasma initial conditions (solid), MC-KLN (dashed), MC-Glauber using participant scaling (dotted) and binary collision scaling (dash-dotted).

As discussed at the outset, MC-KLN fails to describe experimental and simultaneously Esumi (2011); Qiu et al. (2012) because of its small ratio . The fluctuating IP-Glasma initial state presented here has a larger , closer to that of the MC-Glauber model that is tuned to describe experimental reasonably well Schenke et al. (2011c).

In summary, we introduced the IP-Glasma model of fluctuating initial conditions for heavy-ion collisions. This model goes beyond the MC-KLN implementation by using CYM solutions instead of -factorization and including quantum fluctuations on the dynamically generated transverse length scale . Further, unlike MC-KLN, its parameters are fixed by HERA inclusive and diffractive e+p DIS data. At fixed impact parameter, this model naturally produces NBD multiplicity fluctuations that are known to describe and multiplicity distributions, and its ratio of initial triangularity to eccentricity is more compatible with experimental data of harmonic flow coefficients.

Looking forward, an improved matching to the hydrodynamic description, starting at time , can be achieved by including classical Yang-Mills evolution of the system up to this time. However, we do not expect a significant modification of the presented results for and as suggested by previous work Lappi and Venugopalan (2006). Further refinements include treating color charge correlations encoded in the JIMWLK hierarchy for improved rapidity and energy distributions Dumitru et al. (2011b); Iancu and Triantafyllopoulos (2012) and eliminating arbitrariness in choice of thermalization time by an ab initio treatment of thermalization Dusling et al. (2011); Blaizot et al. (2012); Kurkela and Moore (2011); Berges et al. (2012). Detailed studies of higher flow harmonics using dissipative hydrodynamic simulations and comparison to experimental data will allow for further discrimination between different initial conditions. Specifically, it would be interesting to see if these comparisons are able to distinguish between our Glasma flux tube scenario with granularity on the energy dependent scale and other non-perturbative string scenarios which share common features such as NBD fluctuations but are sensitive to  Wang and Gyulassy (1991); Magas et al. (2001); ?; Flensburg et al. (2011).

Acknowledgments We thank S. Chattopadhyay, A. Dumitru, C. Gale, S. Jeon, T. Lappi, L. McLerran, and Z. Qiu for helpful discussions. BPS and RV are supported by US Department of Energy under DOE Contract No.DE-AC02-98CH10886 and acknowledge additional support from a BNL “Lab Directed Research and Development” grant LDRD 10-043.


  1. The IP-Sat model gives good -squared fits to available small HERA data Kowalski et al. (2006) and fixed target e+A DIS data Kowalski et al. (2008). Since the analysis of Ref. Kowalski et al. (2006), more precise data is now available; a repeat analysis is desirable.
  2. Light cone quantities are defined as . The coordinates are defined as and .
  3. The metric in the coordinate system is so that . The components of the gauge field are related by .
  4. As noted previously Blaizot et al. (2010), the CYM approach treats soft modes with more accurately than in commonly used factorized descriptions.
  5. The exact numerical factor between the two quantities depends on the details of the calculation Lappi (2008) but will not be relevant for our final results.
  6. Here, and henceforth, the distributions are evaluated at , for zero rapidity, where is the average transverse momentum of charged hadrons in p+p collisions at a given .
  7. in all computations presented here.
  8. Modifications to Gaussian distributions, relevant for smaller nuclei, have recently been explored in Dumitru et al. (2011a).


  1. B. Alver and G. Roland, Phys. Rev. C81, 054905 (2010).
  2. B. H. Alver, C. Gombeaud, M. Luzum,  and J.-Y. Ollitrault, Phys. Rev. C82, 034913 (2010).
  3. H. Petersen, G.-Y. Qin, S. A. Bass,  and B. Muller, Phys. Rev. C82, 041901 (2010).
  4. B. Schenke, S. Jeon,  and C. Gale, Phys. Rev. Lett. 106, 042301 (2011a).
  5. B. Schenke, S. Jeon,  and C. Gale, Phys. Lett. B702, 59 (2011b).
  6. Z. Qiu and U. W. Heinz, Phys. Rev. C84, 024911 (2011).
  7. Z. Qiu, C. Shen,  and U. W. Heinz, Phys. Lett. B707, 151 (2012).
  8. M. Bleicher, K. Bugaev, P. Rau, A. Sorin, J. Steinheimer, et al.,  (2011), arXiv:1106.3647 [nucl-th] .
  9. B. Schenke, S. Jeon,  and C. Gale, Phys. Rev. C85, 024901 (2011c).
  10. E. Komatsu et al. (WMAP Collaboration), Astrophys.J.Suppl. 180, 330 (2009).
  11. F. Gelis, E. Iancu, J. Jalilian-Marian,  and R. Venugopalan, Ann.Rev.Nucl.Part.Sci. 60, 463 (2010).
  12. A. Kovner, L. D. McLerran,  and H. Weigert, Phys. Rev. D52, 3809 (1995a).
  13. A. Kovner, L. D. McLerran,  and H. Weigert, Phys. Rev. D52, 6231 (1995b).
  14. A. Krasnitz and R. Venugopalan, Phys. Rev. Lett. 84, 4309 (2000).
  15. A. Krasnitz and R. Venugopalan, Phys. Rev. Lett. 86, 1717 (2001).
  16. H. Kowalski, T. Lappi,  and R. Venugopalan, Phys. Rev. Lett. 100, 022303 (2008).
  17. H. I. Miettinen and J. Pumplin, Phys. Rev. D18, 1696 (1978).
  18. S. A. Voloshin, Phys. Lett. B632, 490 (2006).
  19. E. V. Shuryak, Phys. Rev. C76, 047901 (2007).
  20. T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey,  and Y. Nara, Phys. Lett. B636, 299 (2006).
  21. H.-J. Drescher, A. Dumitru, A. Hayashigaki,  and Y. Nara, Phys. Rev. C74, 044905 (2006).
  22. H.-J. Drescher and Y. Nara, Phys. Rev. C76, 041903 (2007).
  23. S. Esumi (PHENIX Collaboration), J.Phys.G G38, 124010 (2011).
  24. J. Bartels, K. J. Golec-Biernat,  and H. Kowalski, Phys. Rev. D66, 014001 (2002).
  25. H. Kowalski and D. Teaney, Phys. Rev. D68, 114005 (2003).
  26. L. D. McLerran and R. Venugopalan, Phys. Rev. D49, 2233 (1994a).
  27. L. D. McLerran and R. Venugopalan, Phys. Rev. D49, 3352 (1994b).
  28. L. D. McLerran and R. Venugopalan, Phys. Rev. D50, 2225 (1994c).
  29. H. Kowalski, L. Motyka,  and G. Watt, Phys. Rev. D74, 074016 (2006).
  30. P. Tribedy and R. Venugopalan, Nucl. Phys. A850, 136 (2011).
  31. P. Tribedy and R. Venugopalan, Phys.Lett. B710, 125 (2012), arXiv:1112.2445 [hep-ph] .
  32. A. Dumitru and Y. Nara,  (2012), arXiv:1201.6382 [nucl-th] .
  33. J. Jalilian-Marian, A. Kovner, L. D. McLerran,  and H. Weigert, Phys. Rev. D55, 5414 (1997).
  34. Y. V. Kovchegov, Phys. Rev. D54, 5463 (1996).
  35. T. Lappi, Phys. Rev. C67, 054903 (2003).
  36. J. P. Blaizot, T. Lappi,  and Y. Mehtar-Tani, Nucl. Phys. A846, 63 (2010).
  37. T. Lappi, Eur. Phys. J. C55, 285 (2008).
  38. A. Krasnitz, Y. Nara,  and R. Venugopalan, Nucl. Phys. A717, 268 (2003).
  39. A. Krasnitz and R. Venugopalan, Nucl. Phys. B557, 237 (1999).
  40. A. Dumitru, J. Jalilian-Marian,  and E. Petreska, Phys. Rev. D84, 014018 (2011a).
  41. P. Romatschke and R. Venugopalan, Phys. Rev. D74, 045011 (2006).
  42. K. Dusling, F. Gelis,  and R. Venugopalan, Nucl.Phys. A872, 161 (2011).
  43. V. Magas, L. Csernai,  and D. Strottman, Phys.Rev. C64, 014901 (2001).
  44. V. Magas, L. Csernai,  and D. Strottman, Nucl.Phys. A712, 167 (2002).
  45. F. Gelis, T. Lappi,  and L. McLerran, Nucl. Phys. A828, 149 (2009).
  46. T. Lappi, S. Srednyak,  and R. Venugopalan, JHEP 01, 066 (2010).
  47. G.-Y. Qin, H. Petersen, S. A. Bass,  and B. Muller, Phys.Rev. C82, 064903 (2010).
  48. T. Lappi and R. Venugopalan, Phys. Rev. C74, 054905 (2006).
  49. T. Lappi, Phys. Lett. B643, 11 (2006).
  50. B. Schenke, S. Jeon,  and C. Gale, Phys. Rev. C82, 014903 (2010).
  51. M. Gyulassy, D. H. Rischke,  and B. Zhang, Nucl.Phys. A613, 397 (1997).
  52. A. Dumitru, J. Jalilian-Marian, T. Lappi, B. Schenke,  and R. Venugopalan, Phys.Lett. B706, 219 (2011b).
  53. E. Iancu and D. Triantafyllopoulos, JHEP 1204, 025 (2012), arXiv:1112.1104 [hep-ph] .
  54. J.-P. Blaizot, F. Gelis, J.-F. Liao, L. McLerran,  and R. Venugopalan, Nucl.Phys. A873, 68 (2012).
  55. A. Kurkela and G. D. Moore, JHEP 1111, 120 (2011).
  56. J. Berges, K. Boguslavski,  and S. Schlichting, Phys.Rev. D85, 076005 (2012), arXiv:1201.3582 [hep-ph] .
  57. X.-N. Wang and M. Gyulassy, Phys.Rev. D44, 3501 (1991).
  58. C. Flensburg, G. Gustafson, L. Lonnblad,  and A. Ster, JHEP 1106, 066 (2011).
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
The feedback must be of minumum 40 characters
Add comment

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question