aa r X i v : . [ h e p - ph ] J u l Glittering Glasmas
F. Gelis a , T. Lappi a,b , L. McLerran a,c a Institut de Physique Th´eorique, Bˆat. 774, CEA/DSM/Saclay, 91191 Gif-sur-Yvette, France b Department of Physics, P.O. Box 35, 40014 University of Jyv¨askyl¨a, Finland c Physics Department and Riken Brookhaven Center, Brookhaven National Laboratory, Upton, NY 11973, USA
Abstract
We compute the production of gluons from Glasma color flux tubes. We calculate the probabilitydistribution of gluon multiplicities arising from the distribution of color electric and color magneticflux tubes found in the Glasma. We show that the result corresponds to the negative binomialprobability distribution observed in experiments. The parameter k that characterizes this distri-bution is proportional to the number of colors N c2 − Key words: glasma, multiplicity distribution
PACS:
1. Introduction
In high energy nuclear collisions, it has been argued that a Glasma is formed through thecollision of two sheets of Color Glass Condensate (CGC) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,14, 15, 16, 17, 18] (for a review and additional references see [19, 20]). At very high energies,the gluon density per unit area per unit rapidity, d N/ d y d r ⊥ , is of order Q /α s in the CGC andGlasma. This density is so large that the interaction strength of QCD is weak, α s ≪
1. Indeed,in these relations, the scale of the coupling constant is set by the saturation momentum Q s , andthe saturation momentum grows with both increasing energy and size of the nucleus. The Glasmais formed during the time it takes two Lorentz contracted sheets of Colored Glass to pass throughone another, τ ∼ e − κ/α s /Q s a time parametrically short compared to the natural time scale fordecay of the flux tubes, τ decay ∼ /Q s . We shall not describe in detail the properties of either theGlasma or the CGC in this paper and refer the reader to the original literature for details.In this paper, we compute the probability distribution for the multiplicity of gluons producedin the Glasma. We show that the distribution of gluons arising from such a decay is not a Poissondistribution as might be expected for the decay of an external source into particles in a weaklyinteracting theory. It turns out the distribution is a negative binomial distribution.Recall that a Poisson distribution, P Poisson n = 1 n ! ¯ n n e − ¯ n (1)is completely characterized by its mean value ¯ n . In contrast, a negative binomial distribution is a2-parameter distribution of the form P NB n = Γ( k + n )Γ( k )Γ( n + 1) ¯ n n k k (¯ n + k ) n + k . (2) Preprint submitted to Elsevier his distribution has larger fluctuations than a Poisson distribution, but tends to a Poisson distri-bution if k → + ∞ at fixed ¯ n .For the processes we consider, the decays will be into one coherent state associated with agluon. For this reason, and because a negative binomial distribution does not fall as 1 /n ! at large n , as does a Poisson distribution, we will refer to this property of the distribution as tenacious. Thetenaciousness of this distribution results in an amplification of the intensity of multiply emittedgluons relative to that of a Poisson distribution. This means that there are larger fluctuations inhigh gluon multiplicity events than would be typical of a Poisson distribution. Thus, we shall usethe acronym “Glitter” to describe radiation from these flux tubes, as an abbreviation for GLuonIntensification Through Tenacious Emission of Radiation.In this paper, we shall compute the multiplicity distribution of gluons produced by Glasma fluxtubes. We show that a single flux tube decays into gluons with a negative binomial distributionand is characterized by a parameter k of order one. This implies that a decay of N FT flux tubesproduces a negative binomial distribution characterized by k = N FT k . For k = 1, a negativebinomial distribution is a Bose-Einstein distribution , so single flux tube decays are close in formto Bose-Einstein distributions. In addition to the parameter k , the negative binomial distributionis also parameterized by the average multiplicity ¯ n . The ratio ¯ n/k = ¯ n/N FT k is approximatelythe multiplicity per flux tube.Based on this interpretation of the decay of an ensemble of flux tubes we argue that the fluxtubes are a “glittering” glasma. We will then iscuss our results in the context of the extraction ofnegative binomial distribution from the UA(5) and PHENIX collaborations.
2. Calculation of the multiplicity distribution
The probability distribution of a discrete quantity is conveniently defined in terms of its gen-erating function F ( z ) ≡ ∞ X n =0 z n P n , (3)from which one can compute the moments of the multiplicity distribution as h n ( n − · · · ( n − q + 1) i = d q F ( z )d z q (cid:12)(cid:12)(cid:12)(cid:12) z =1 . (4)(The moments defined in this way are known as the factorial moments.) It was shown in Refs. [21,22, 23] that in the case of the central rapidity region of nucleus-nucleus collisions, when both nucleican be described as strong color sources ρ ∼ /g , these moments can be computed as h n ( n − · · · ( n − q + 1) i = Z [d ρ ][d ρ ] W y [ ρ ] W y [ ρ ] (cid:0) n [ ρ , ρ ] (cid:1) q . (5)This result is valid to leading log accuracy, i.e. it includes the leading order in α s with all thepowers of α s ln 1 /x resummed into the rapidity dependence of the weight functionals W y [ ρ ]. The A Bose-Einstein, or geometrical, distribution is a thermal distribution for single state systems. Its probabilitydistribution reads: P BE n = 11 + ¯ n „ ¯ n n « n . p-kk = p Figure 1: The fundamental building block, the one gluon production amplitude. p p ... p q p q ... p p Figure 2: Diagrams that have to be contracted into connected one. factor n [ ρ , ρ ] inside the integral in the r.h.s. of eq. (5) is the integrated multiplicity correspondingto a fixed configuration of color charge densities ρ and ρ : n [ ρ , ρ ] = Z d p ⊥ d y p d N d p ⊥ d y p , (6)where d N / d p ⊥ d y p is obtained by Fourier transforming the classical gauge field radiated by thesources ρ , . In our leading log calculation the classical fields and thus also the single gluonmultiplicity are boost invariant, so the integral over rapidity in eq. (6) is just a constant factor.We are assuming that the rapidity interval is small enough compared to 1 /α s ; otherwise there areadditional large logarithms that must be resummed; this case is studied in detail in Ref. [23].In the Glasma the single inclusive multiplicity is of order 1 /α s and thus the moment defined ineq. (5) is of order (1 /α s ) q . We are computing the moments only to leading order in α s , thus powersof n lower than q are negligible compared to the q ’th moment. In particular h n ( n − · · · ( n − q +1) i ≈h n q i at this level of accuracy. This means that the leading (in α s ) correlations in the multiplicitydistribution come entirely from the average over the distribution of sources ρ , ; i.e. the dominantbehavior of the probability distribution comes from the large logarithms of the energy resummedinto the W ’s. The contributions to the probability distribution for fixed sources from higher looporders that were studied in Ref. [24] are suppressed by powers of α s and thus contribute to ourcalculation only when they are enhanced by large logarithms.The fundamental properties of the probability distribution are better reflected in the factorialcumulants m q ≡ h n ( n − · · · ( n − q + 1) i − disc. = d q ln F ( z )d z q (cid:12)(cid:12)(cid:12)(cid:12) z =1 , (7)where “disc.” denotes the disconnected contributions that can be expressed in terms of the lowercumulants. We shall now turn to calculating these quantities of the multiplicity distribution in the Note that the factorial cumulant defined here differs from the conventional one in that we are defining thecumulant from h n ( n − · · · ( n − q + 1) i instead of h n q i . In terms of the generating function this is a question ofdifferentiating w.r.t. z in stead of ln z . As explained above the difference between the two is of higher order in α s than our calculation and we will not discuss it further. -k k -k kp-k k-p Figure 3: Contraction contributing to the dominant correlation, building block of rainbow diagram p p p -k k-p p -kk-p k -k Figure 4: Contraction contributing to a subdominant correlation, non-rainbow diagram
Glasma at the lowest nontrivial order in the sources. This approximation is equivalent to assumingthat the momenta of the produced gluons are all larger than the saturation scale. For concretenessone can take the distribution of the sources from the MV model: W [ ρ ] = C exp (cid:20) − Z d x ⊥ ρ a ( x ⊥ ) ρ a ( x ⊥ ) g µ (cid:21) . (8)Since the probability distribution essentially depends only on the combinatorics of pairwise sourceconnections, our result applies equally well to a nonlocal Gaussian distribution that would moreclosely reproduce a solution of the BK equation [25, 4, 26].The calculation proceeds in the same way as that of the second and third cumulants computedin Refs.[27, 28, 29], and we refer to these works for a more detailed description. Computing theprobability distribution to all orders in the sources ρ , is in principle possible using the methodsdeveloped for the single inclusive gluon production [14, 15, 16, 17], but reliably computing thehigher cumulants requires a significant numerical effort to gather enough statistics and is left forfuture work.The fundamental building block in our calculation is the amplitude to produce one gluon withmomentum p from the fixed classical color charges ρ ( k ⊥ ) and ρ ( p ⊥ − k ⊥ ). This amplitude reads ρ ( k ⊥ ) k ⊥ ρ ( p ⊥ − k ⊥ )( p ⊥ − k ⊥ ) L γ ( p , k ⊥ ) , (9)where we are not writing the color indices explicitly. Here, L γ ( p , k ⊥ ) denotes the effective Lipatovvertex. We do not need the explicit expression of its components, and it will be sufficient to knowthat it satisfies the following two properties L γ ( p , k ⊥ ) = L γ ( p , p ⊥ − k ⊥ ) ,L γ ( p , k ⊥ ) L γ ( p , k ⊥ ) = − p ⊥ − k ⊥ ) k ⊥ p ⊥ . (10)4 ... p q p q ... p Figure 5: Rainbow diagram. p p = p Figure 6: Dimer notation for rainbow-like links.
The diagrammatic notation for this amplitude is shown in fig. 1. To compute the q ’th cumulantwe need to take 2 q factors of this basic building block ( q for the amplitude and q for the complexconjugate) and perform the averages over the sources. Because the distribution of the sources ineq. (8) is Gaussian, we only need to keep track of contractions of pairs of sources ρ and (separately)of pairs of sources ρ and replace them by the correlator h ρ ( k ⊥ ) ρ ( k ′⊥ ) i = (2 π ) δ ( k ⊥ + k ′⊥ ) g µ ( k ⊥ ) . (11)For the time being we shall leave an unspecified k ⊥ -dependence in the correlation function g µ ( k ⊥ ).In the MV model [1, 2, 3] g µ ( k ⊥ ) is a constant, but JIMWLK or BK evolution can effectivelylead to a different k ⊥ -dependence [30]. With the simplified diagrammatic notation introduced infig. 1 this combinatoric problem now corresponds to forming a connected contraction of the 2 q boxes, each with two lines attached, illustrated in fig. 2.We shall now show why the only contributing contractions are “rainbow” diagrams, where onone side (upper or lower) of the diagram the two boxes corresponding to the same momentum p r are contracted with each other, as in fig. 5. An example of a building block of a rainbow diagramis given in fig. 3. There are a total of four propagators with the same momentum k ⊥ . Two ofthese are cancelled by the contractions of the four Lipatov vertices attached to the ends, leavinga quadratically infrared divergent contribution ∼ d k ⊥ / k ⊥ to the integral over k ⊥ . These kindsof divergences are a sign that this contribution to the multiparticle correlation is sensitive to thewhole correlated area in the transverse plane. They are regulated at the scale Q s (since Q − is thecorrelation length between color charges in the transverse direction), giving a contribution of order1 /Q . Compare this to the “non-rainbow” contribution depicted in fig. 4. Here there are only twopropagators with the same momenta and the k ⊥ -integral is convergent. Instead of O (1 /Q ), thisyields a contribution ∼ / p ⊥ or ∼ / p ⊥ which shall be neglected here since we are assuming p ⊥ ≫ Q s . Thus the contractions of the boxes in fig. 2 have to form a “rainbow diagram” (fig. 5)in either the upper or lower part of the diagram. There cannot be a rainbow on both sides sincethis would lead to a disconnected contribution.Now that we have reduced the combinatoric problem to rainbow diagrams we can simplify our5 p p p Figure 7: Connected diagram in the polymer notation introduced in fig. 6 diagrammatic notation even further, as shown in fig. 6. We can consider the two boxes connected bythe line in the rainbow as a dimer with two ends. Now we must count the number of ways in whichour q dimers can be combined to form a connected loop. The first end of the first dimer can beconnected to 2 q − r . The second endof dimer r now has 2 q − q − · (2 q − ·· · ·· q − ( q − q ( q − q = 2 and q = 3 this gives4 and 16, which were the number of diagrams evaluated in [27] and [29].The color structure can be simplified in the following way. Let us denote the color indices ofthe sources i = 1 . . . q contracted on the upper side of the diagram in fig. 5 by a i (because of therainbow structure of the connections on the upper side, we do not need separate indices for thesources in the complex conjugate amplitude), the color index of the gluon with p i by b i and thecolor indices on the lower side of the diagram c i (source connected to gluon p i in the amplitude)and c ′ i (in the complex conjugate). Now the color structure of the upper side of the rainbow and theboxes is f a b c f a b c ′ . . . f a q b q c q f a q b q c ′ q = ( C A ) q δ c c ′ . . . δ c q c ′ q . C A = N c is the Casimir of the adjointrepresentation. The c i , c ′ i indices must now be contracted pairwise into a single connected loop asin fig. 7, which yields a factor tr(1 adj ) = N c2 −
1, making the total color factor ( N c ) q ( N c2 − q transverse momentum integrals from the powers of the sources.There are 2 q delta functions from the source correlators eq. (11) and 2 q momentum conservationdelta functions from the three gluon vertices. Not all these delta functions are independent: twoof them end up having the same argument, yielding one factor of the transverse area (denoted S ⊥ ). Therefore, there is only one remaining transverse momentum integral. One can choose thisremaining momentum to be the one circulating in all of the lower part of the diagram, which weshall denote by k ⊥ (this is the momentum that circulates along the loop in fig. 7). On the rainbowside of the diagram there is a squared propagator 1 / ( k ⊥ − p ⊥ i ) for all the sources i = 1 . . . q . Onthe non-rainbow side the transverse momentum in all the propagators is the same, giving a factor1 / k q ⊥ . Half of these propagators are cancelled by the squares of the Lipatov vertices, which alsocontribute an inverse square of the external momentum. Combining the combinatorial factors fromthe source averages, the propagators, Lipatov vertices and factors from the invariant measure, we6nally get * d N d y d p ⊥ . . . d y q d p ⊥ q + conn. = h q ( q − i ( N c ) q ( N c2 − S ⊥ ( p ⊥ ) · · · ( p ⊥ q ) g q q (2 π ) q × Z d k ⊥ (2 π ) (cid:18) g µ ( k ⊥ ) k ⊥ (cid:19) q g µ ( p ⊥ − k ⊥ )( p ⊥ − k ⊥ ) · · · g µ ( p ⊥ q − k ⊥ )( p ⊥ q − k ⊥ ) . (12)This general formula also reproduces the result of refs. [12, 31, 32] for the single inclusive spectrumcase q = 1; in this case the combinatorial factor in the square bracket must be taken to be 1 insteadof 2 to avoid double counting the only contributing diagram.The weak source result eq. (12) is infrared divergent in the MV model ( g µ ( k ⊥ ) constant).Physically this is modified by several effects. Even in the weak field limit BK or BFKL evolutionleads to an anomalous dimension 0 < γ < g µ ( k ⊥ ) ∼ k − γ ) ⊥ in the geometric scaling region k ⊥ & Q s . Deep in the saturation regime it has been argued [30]that the correlator effectively behaves as g µ ( k ⊥ ) ∼ k ⊥ . Ultimately the infrared behavior of themultigluon spectrum is regulated by the nonlinear interactions that are not included in our presentcomputation. This is seen explicitly and analytically in the “pA” case [33, 34, 35] and in numericalcomputations of the glasma fields in the fully nonlinear case [14, 15, 16, 17]. Since the full nonlineardynamics are known to regulate the infrared behavior in the case of the single gluon spectrum wehave strong reasons to expect that they will also do so in the case of multiple gluon production;at the same scale k ⊥ . Q s . We emphasize that an essential point in this argument is that thequantity appearing in eq. (12) is not a single color charge correlator divided by a large power k q ⊥ ,but the same correlator g µ ( k ⊥ ) / k ⊥ that appears in the single inclusive gluon spectrum raised toa large power q .The effect of saturation on the multigluon spectrum at k ⊥ . Q s has a very intuitive interpreta-tion in the the glasma flux tube picture. The size of the flux tube, 1 /Q s , is the correlation length ofthe system and we should not have contributions from longer distance scales. We effectively takethis into account by regulating all the infrared divergences at the scale Q s , and thus approximatingthe integral in eq. (12) by (cid:16) πκ q − Q q − p ⊥ · · · p ⊥ q (cid:17) − . Here κ is a constant of order one thatdepends on the details of how the infrared divergences are regulated at the scale Q s . Since Q s isthe typical momentum of the produced gluons, not a lower limit, we expect that numerically κ < : (cid:28) d N d y d p ⊥ (cid:29) ≈ N c ( N c2 − π g S ⊥ ( g µ ) p ⊥ . (13) In comparing this to the result in [27] (eq. (27)) note that there is an erroneous factor of 1 / − q . In addition, there is an overall error of (2 π ) . These errors werepresent also in Ref. [29], but have been corrected in the published version. We are neglecting an additional logarithmic dependence in p ⊥ .
7e can now express our result as * d N d y d p ⊥ . . . d y q d p ⊥ q + conn. = ( q − N c2 − κQ S ⊥ π (cid:18) ( g µ ) g π N c κQ (cid:19) q p ⊥ ) · · · ( p ⊥ q ) = ( q − N c2 − κQ S ⊥ π D d N d y d p ⊥ E . . . D d N d y q d p ⊥ q E(cid:0) ( N c2 − κQ S ⊥ / (2 π ) (cid:1) q (14)If we integrate this equation over the rapidities and transverse momenta of the q gluons, againconsistently regulating all the infrared divergences at the scale Q s , we obtain our result for thefactorial cumulant as m q = ( q − k (cid:16) ¯ nk (cid:17) q , (15)with k = κ ( N c2 − Q S ⊥ π . (16)The exact constant factors, encoded in the coefficient κ , depend on the exact way the infrareddivergences (logarithmic for the single inclusive, power law for the multigluon correlations) areregulated. These factors cannot be obtained exactly in an analytic calculation to the lowest orderin the sources. However, the main parametric dependences in the relevant variables α s , Q s , S ⊥ and N c can be expected to be the same to all orders in the sources. A possible additional (mild) q -dependence in κ would be a minor correction to the behavior of the probability distribution,mostly determined by the combinatorial factor ( q − negative binomial distribution. It arises very naturallyin the Glasma based on the Gaussian combinatorics of the classical sources and the assumption ofthe fluctuations in the system being dominated by a correlation length 1 /Q s .
3. Glittering Glasma: Interpretation of the result
A negative binomial distribution is characterized by two parameters, the mean ¯ n and k , interms of which the probability to produce n particles is P NB n = Γ( k + n )Γ( k )Γ( n + 1) ¯ n n k k (¯ n + k ) n + k . (17)The distribution is characterized by the generating function F k, ¯ n ( z ) ≡ ∞ X n =0 z n P n = (cid:16) − ¯ nk ( z − (cid:17) − k . (18)The moments of the distribution can be obtained from the generating function by differentiatingwith respect to z . The connected parts of the moments, or cumulants, are generated by the loga-rithm of the generating function. The factorial cumulants (7) of the negative binomial distributionare given by m q ≡ d q d z q ln F k, ¯ n ( z ) (cid:12)(cid:12)(cid:12)(cid:12) z =1 = ( q − k (cid:16) ¯ nk (cid:17) q . (19)8n contrast, for a Poisson distribution m = ¯ n and m q = 0 for q >
1. This quantity is theexpectation value m q = h n ( n − · · · ( n − q + 1) i − disc. , (20)where the disconnected part “disc” can be expressed in terms of the lower order cumulants. Theconvention that the expectation value in (20) is taken of the product n ( n − · · · ( n − p + 1) andnot of n p means that we are subtracting the “Poissonian” part from the moment, which is whyall the factorial cumulants of a Poisson distribution are zero for q ≥
2. Note that the Poissonianpart is suppressed by powers of α s in the Glasma; thus in practice the difference between m p anda conventional cumulant is neglected in our analysis.Two common special cases of the negative binomial are the Poisson distribution, obtained inthe limit k → ∞ at fixed ¯ n , and the geometrical–or Bose-Einstein–distribution obtained when k = 1. The negative binomial distribution is wider than a Poisson distribution typically associatedwith independent emission of particles; this can be seen e.g. from the variance σ = h n i − h n i = ¯ n + ¯ n k . (21)A useful property of the negative binomial distribution with parameters ¯ n, k is that it is is alsothe distribution of a sum of k independent random variables drawn from a Bose-Einstein distribu-tion with mean ¯ n/k . This is easily seen from the generating function in eq. (18), remembering thatthe generating function of a sum of independent random variables is the product of their generat-ing functions . This has a consequence that an incoherent superposition of N emitters that havea negative binomial distribution with parameters k , ¯ n produces a negative binomial distributionwith parameters N k , N ¯ n .A natural physical interpretation of our result can be given in terms of emission from indepen-dent glasma flux tubes. Geometrically, the transverse area S ⊥ is filled with Q S ⊥ independent fluxtubes of size ∼ /Q . Each of these tubes emits gluons in N c2 − k increases somewhat with δη , the size of the rapidity interval in which the particles are measured.The dependence is, however, very slow for large δη , pointing to the presence of a long rangecorrelation in the system [43]. This is natural in the Glasma picture, since flux tubes extend overlarge rapidity intervals. The number of flux tubes, which gives the parameter k of the negative Explicitly, consider n = n + · · · + n r where the n i ’s are independent of each other. The probability dis-tribution of n is then P n = P n · · · P n r δ ( n − P ri =1 n i ) P n · · · P n r and the generating function P n z n P n = P n · · · P n r z n + ··· + n r P n · · · P n r , which is the product of the individual generating functions for the variables n i . p ¯ p , thepicture we present seems fairly consistent. For a fixed collision energy we would expect scaling k ∼ Q S ⊥ ∼ N part . While keeping this caveat in mind, the results from UA5 [39] and E735 [45]( k ≈ . . . N part = 2) and PHENIX k ≈
350 for 0-5% most central ( N part ≈ k = 690 when extrapolated to a zero centrality bin width [42], seem very consistent with thisestimate.From the experimental fit of the parameter k in central gold-gold collisions by PHENIX [42]and the value Q s ≈ . κ ≈ . κ . This, however, means that at RHIC energies the flux tube size, as measured in themultiplicity distribution, is not yet very clearly separated from the confinement scale. At LHCenergies we can expect this separation to be clearer.For increasing collision energy we would expect Q s and therefore k to increase. The energieswhere the UA5 measurements are done are still in the transition region from a behavior of k decreasing with energy from lower √ s , but we would expect k at the LHC to be clearly larger.This decreasing behavior at low energy follows because of the Poisson nature of low energy particleemission, and that for a Poisson distribution k → ∞ .The negative binomial has been interpreted as resulting from a partial stimulated emissionor cascade process [43]. It has been known in the literature [46, 47] that the distribution wouldnaturally arise from a superposition of subsystems with Bose-Einstein distributions. Nevertheless,a popular approach has remained to interpret the observations in terms of a fluctuating numberof strings [48], each producing particles typically with a Poisson distribution [49, 50] (see also[51, 52, 53] for a more pQCD based approach). While the picture of flux tubes in the glasma hasmany similarities to ideas in string model phenomenology, the distribution of particles producedfrom one flux tube is different. The probability distribution of gluons from a glasma flux tube isnot a narrow Poissonian, but has very large fluctuations: the glittering of the glasma.
4. Summary
The Glasma provides a successful phenomenology of a particle production in high energyhadronic collisions. There is now experimental data on the ridge phenomena that show fluxtube structures in two particle correlations [54, 55, 56]. In addition, long range correlations ofremarkable strength are seen in heavy ion collisions [44].The Glitter of the flux tube decay may provide a strong tool for disentangling various de-scriptions of the flux tubes, since it naturally leads to a negative binomial distribution for themultiplicity of produced particles. However, in order to make a more convincing case for the originof the negative binomial distribution of particle multiplicities, we need a systematic study of binsize effects on the extraction of the parameter k from various centralities of heavy ion collisions.10 cknowledgements The authors gratefully acknowledge conversations with Raju Venugopalan. L. McLerran wassupported in part by the Theoretical Physics Division at CEA-Saclay, and this work is a product ofthe stimulating intellectual atmosphere there. The research of L. McLerran is supported under DOEContract No. DE-AC02-98CH10886. T. Lappi is supported by the Academy of Finland, project126604. F. Gelis is supported in part by Agence Nationale de la Recherche via the programmeANR-06-BLAN-0285-01.
References [1] L. D. McLerran and R. Venugopalan, Phys. Rev.
D49 , 2233 (1994), [arXiv:hep-ph/9309289].[2] L. D. McLerran and R. Venugopalan, Phys. Rev.
D49 , 3352 (1994), [arXiv:hep-ph/9311205].[3] L. D. McLerran and R. Venugopalan, Phys. Rev.
D50 , 2225 (1994), [arXiv:hep-ph/9402335].[4] Y. V. Kovchegov, Phys. Rev.
D54 , 5463 (1996), [arXiv:hep-ph/9605446].[5] J. Jalilian-Marian, A. Kovner, L. D. McLerran and H. Weigert, Phys. Rev.
D55 , 5414 (1997),[arXiv:hep-ph/9606337].[6] J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, Nucl. Phys.
B504 , 415 (1997),[arXiv:hep-ph/9701284].[7] J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, Phys. Rev.
D59 , 014014 (1999),[arXiv:hep-ph/9706377].[8] E. Iancu, A. Leonidov and L. D. McLerran, Nucl. Phys.
A692 , 583 (2001), [arXiv:hep-ph/0011241].[9] E. Ferreiro, E. Iancu, A. Leonidov and L. McLerran, Nucl. Phys.
A703 , 489 (2002), [arXiv:hep-ph/0109115].[10] E. Iancu and L. D. McLerran, Phys. Lett.
B510 , 145 (2001), [arXiv:hep-ph/0103032].[11] A. H. Mueller, Phys. Lett.
B523 , 243 (2001), [arXiv:hep-ph/0110169].[12] A. Kovner, L. D. McLerran and H. Weigert, Phys. Rev.
D52 , 3809 (1995), [arXiv:hep-ph/9505320].[13] A. Kovner, L. D. McLerran and H. Weigert, Phys. Rev.
D52 , 6231 (1995), [arXiv:hep-ph/9502289].[14] A. Krasnitz and R. Venugopalan, Nucl. Phys.
B557 , 237 (1999), [arXiv:hep-ph/9809433].[15] A. Krasnitz and R. Venugopalan, Phys. Rev. Lett. , 1717 (2001), [arXiv:hep-ph/0007108].[16] A. Krasnitz, Y. Nara and R. Venugopalan, Phys. Rev. Lett. , 192302 (2001), [arXiv:hep-ph/0108092].[17] T. Lappi, Phys. Rev. C67 , 054903 (2003), [arXiv:hep-ph/0303076].[18] T. Lappi and L. McLerran, Nucl. Phys.
A772 , 200 (2006), [arXiv:hep-ph/0602189].[19] E. Iancu and R. Venugopalan, The color glass condensate and high energy scattering in QCD, in
Quark gluonplasma , edited by R. Hwa and X. N. Wang, World Scientific, 2003, arXiv:hep-ph/0303204.[20] H. Weigert, Prog. Part. Nucl. Phys. , 461 (2005), [arXiv:hep-ph/0501087].[21] F. Gelis, T. Lappi and R. Venugopalan, Phys. Rev. D78 , 054019 (2008), [arXiv:0804.2630 [hep-ph]].[22] F. Gelis, T. Lappi and R. Venugopalan, Phys. Rev.
D78 , 054020 (2008), [arXiv:0807.1306 [hep-ph]].[23] F. Gelis, T. Lappi and R. Venugopalan, Phys. Rev.
D79 , 094017 (2008), [arXiv:0810.4829 [hep-ph]].[24] F. Gelis and R. Venugopalan, Nucl. Phys.
A776 , 135 (2006), [arXiv:hep-ph/0601209].[25] I. Balitsky, Nucl. Phys.
B463 , 99 (1996), [arXiv:hep-ph/9509348].[26] Y. V. Kovchegov, Phys. Rev.
D61 , 074018 (2000), [arXiv:hep-ph/9905214].[27] A. Dumitru, F. Gelis, L. McLerran and R. Venugopalan, Nucl. Phys.
A810 , 91 (2008),[arXiv:0804.3858 [hep-ph]].[28] S. Gavin, L. McLerran and G. Moschelli, Phys. Rev.
C79 , 051902 (2009), [arXiv:0806.4718 [nucl-th]].[29] K. Dusling, D. Fernandez-Fraile and R. Venugopalan, arXiv:0902.4435 [nucl-th].[30] E. Iancu, K. Itakura and L. McLerran, Nucl. Phys.
A724 , 181 (2003), [arXiv:hep-ph/0212123].[31] Y. V. Kovchegov and D. H. Rischke, Phys. Rev.
C56 , 1084 (1997), [arXiv:hep-ph/9704201].[32] M. Gyulassy and L. D. McLerran, Phys. Rev.
C56 , 2219 (1997), [arXiv:nucl-th/9704034].[33] A. Dumitru and L. D. McLerran, Nucl. Phys.
A700 , 492 (2002), [arXiv:hep-ph/0105268].[34] D. Kharzeev, Y. V. Kovchegov and K. Tuchin, Phys. Rev.
D68 , 094013 (2003), [arXiv:hep-ph/0307037].[35] J. P. Blaizot, F. Gelis and R. Venugopalan, Nucl. Phys.
A743 , 13 (2004), [arXiv:hep-ph/0402256].[36] UA1, G. Arnison et al. , Phys. Lett.
B123 , 108 (1983).[37] UA5, G. J. Alner et al. , Phys. Lett.
B160 , 193 (1985).[38] UA5, G. J. Alner et al. , Phys. Lett.
B160 , 199 (1985).
39] UA5, R. E. Ansorge et al. , Z. Phys.
C37 , 191 (1988).[40] A. Giovannini and R. Ugoccioni, Int. J. Mod. Phys.
A20 , 3897 (2005), [arXiv:hep-ph/0405251].[41] PHENIX, S. S. Adler et al. , Phys. Rev.
C76 , 034903 (2007), [arXiv:0704.2894 [nucl-ex]].[42] PHENIX, A. Adare et al. , Phys. Rev.
C78 , 044902 (2008), [arXiv:0805.1521 [nucl-ex]].[43] A. Giovannini and L. Van Hove, Z. Phys.
C30 , 391 (1986).[44] STAR, B. K. Srivastava, Int. J. Mod. Phys.
E16 , 3371 (2007), [arXiv:nucl-ex/0702054].[45] E735, C. S. Lindsey et al. , Nucl. Phys.
A544 , 343 (1992).[46] M. Biyajima, N. Suzuki, G. Wilk and Z. Wlodarczyk, Phys. Lett.
B386 , 297 (1996), [arXiv:hep-ph/9507210].[47] G. Wilk and Z. Wlodarczyk, Physica
A376 , 279 (2007), [arXiv:cond-mat/0603157].[48] B. Andersson, G. Gustafson, G. Ingelman and T. Sjostrand, Phys. Rept. , 31 (1983).[49] M. A. Braun, C. Pajares and V. V. Vechernin, Phys. Lett. B493 , 54 (2000), [arXiv:hep-ph/0007241].[50] J. Dias de Deus, E. G. Ferreiro, C. Pajares and R. Ugoccioni, Eur. Phys. J.
C40 , 229 (2005),[arXiv:hep-ph/0304068].[51] I. M. Dremin, Phys. Lett.
B313 , 209 (1993).[52] I. M. Dremin and J. W. Gary, Phys. Rept. , 301 (2001), [arXiv:hep-ph/0004215].[53] I. M. Dremin and V. A. Nechitailo, Phys. Rev.
D70 , 034005 (2004), [arXiv:hep-ph/0402286].[54] J. Putschke, J. Phys.
G34 , S679 (2007), [arXiv:nucl-ex/0701074].[55] STAR, M. Daugherity, J. Phys.
G35 , 104090 (2008), [arXiv:0806.2121 [nucl-ex]].[56] PHOBOS, B. Alver et al. , J. Phys.
G35 , 104080 (2008), [arXiv:0804.3038 [nucl-ex]]., 104080 (2008), [arXiv:0804.3038 [nucl-ex]].