Consistent Determination of Particle Production in p-A, d-A and A-A Collisions using Color Glass Condensate
The success of nonviscous hydrodynamics in describing the collective flow properties of bulk low observables at RHIC has led to the claim that a novel form of strongly coupled Quark Gluon Plasma (sQGP) is created in 200 AGeV Au+Au collisions. This success depends strongly, however, on the initial conditions assumed in the calculations. In particular, agreement with data is only obtained assuming Glauber nuclear reaction plane (participant) geometry. The KLN model of Color Glass Condenstate (CGC) initial conditions require the existence of nonvanishing viscous effects. We develop an improved model of CGC that is more internally consistent between the central and diffuse edge regions. The improved model is shown to describe bulk rapidity distributions for a wide assortment of system sizes, geometries and energies. The self consistency forces a specific dependence of the multiplicity which leads to surprisingly low particle production predictions for Pb+Pb collisions at LHC COM energies. These predictions are similar to those made using simple linear extrapolations of multiplicity attempted by the PHOBOS collaboration.
pacs:12.38.Mh; 24.85.+p; 25.75.-q
One of the most exciting discoveries in high energy nucleus-nucleus collisions at the Relativistic Heavy Ion Collider (RHIC) is the large magnitude of the elliptic flow parameter . characterizes the transverse azimuthal distribution . The impact parameter or participant nucleon dependence of provides a clear signature of the strong collective behavior of the created medium (1); (2); (3); (4). The systematics of (with energy, nuclear size, and centrality) measured at RHIC is consistent with a description of the created Quark Gluon Plasma (QGP) core (with MeV) as a non-viscous “perfect” fluid that evolves under Euler hydrodynamic equations of motion from an early thermalization time fm (5); (6) until the hadronization hypersurface () is reached. Below the highly viscous hadronic resonance gas “corona” evolves with conventional hadronic transport (7); (8). Such a hybrid perfect fluid sQGP core plus nonequilibrium hadron resonance gas corona description, however, is particularly sensitive to the assumed initial geometric distribution of partons produced in nuclear collisions since . Here, the averaging is over the transverse plane where and are the transverse spatial coordinates and the initial number density is used as the weight. The success of the hybrid approach to describe the data is contingent on the validity of the widely assumed Glauber model of participant nucleon geometry.
With Glauber initial conditions, the parton interaction cross sections needed to generate the collective behavior seen at RHIC are much larger than can be supported via purely perturbative collisional processes (9). This has led to the recent strongly coupled (s)QGP paradigm, which is supported by the success of gauge-string dual Anti-deSitter/Conformal Field Theory (AdS/CFT) (10) models in describing the low shear viscosity implied by RHIC data (11); (12); (13).
Hydrodynamic equations of motion, of course, must be supplemented by specified initial temperature and flow velocity fields. Conventionally, these are assumed to be specified on a thermalization proper time fm hypersurface (5); (6). The agreement with data mentioned previously is achieved using the optical Glauber model to calculate the initial distribution of produced matter in the plane transverse to the beam direction in conjunction with Brodsky-Gunion-Kuhn (BGK) (14); (15) dynamics that distribute initial matter in a peculiar “trapezoidal” manner in the longitudinal rapidity direction (8); (16).
There exist, however, alternative and theoretically preferable models for predicting the initial conditions of ultrarelativistic heavy ion reactions that take into account the physics of gluon saturation at small Bjorken . Color Glass Condensate (CGC) refers to a set of models that aim to calculate QCD interactions at asymptotically large energies (at very small Bjorken ) where gluonic degrees of freedom are dominant and the large occupation numbers of gluons provide a classical field description of the physics; the so-called gluon saturation regime (17); (18); (19); (20). Calculations of physical interest (e.g. unintegrated gluon distribution functions, total cross sections) can be attempted perturbatively because the characteristic scale at which non linearities become important, , provides a large momentum scale relative to .
An implementation of the Kharzeev-Levin-Nardi (KLN) (21) parameterization of the unintegrated gluon distributions from CGC is used in (7); (8) to show that CGC/KLN initial conditions, in conjunction with ideal hydrodynamics in the QGP stage, severely overpredict the low data at RHIC. This is explained by the fact that the CGC approach generically predicts a significantly larger initial spatial eccentricity of the produced matter than Glauber/BGK models (22); (23). This larger initial spatial eccentricity is efficiently transferred by ideal hydrodynamic motion into a large elliptical anisotropy in momentum space; a large . Thus, if the CGC initial conditions are the ones that actually reign in HIC at RHIC, dissipative effects would have to be included in QGP evolution because hadronic dissipation has already been taken into account (8). Perturbative QCD cross sections could then possibly generate the amount of collectivity required at RHIC.
There remain, however, serious conceptual problems in the original KLN implementation of the CGC. As pointed out in (22); (23); (26); (24), the definition of in the model is non-universal because it is defined to be proportional to the Glauber participant density, which implicitly depends on the partner nucleus. This same problem leads to a drop in the nuclear saturation momentum swiftly to zero as one goes to the nuclear edge, as opposed to the correct lower limit of a single nucleon saturation momentum (22); (23). There are also further problems stemming from the uncertainty in the evolution and proper treatment of the saturation momentum as one increases Bjorken (26). The detailed properties of the local as one progresses to the nuclear edge are crucial to the determination of the correct initial state predictions from CGC, especially in peripheral collisions.
In this paper, we extend the factorized KLN (fKLN) approach introduced in (23) (what we call the extended local fKLN (elfKLN)) in order to correctly calculate the local density evolution of and consistently calculate the multiplicity of charged particles, , for multiple systems and multiple collision energies. We take special care to account for a consistent evaluation of the model for different system types and sizes and show that this reproduces the multiplicities in both A-A and d-A collisions at RHIC with no varying of input parameters. The self consistency conditions force a particular COM energy dependence on the multiplicities calculated from the model which lead to surprisingly low predictions for Pb+Pb collisions at LHC energies. These predictions are similar to ones made from simple linear extrapolations with energy of particle production from the PHOBOS collaboration (27). We do not consider the high Bjorken evolution of as it does not effect determinations of bulk multiplicity and is better probed with more differential obervables, e.g. detailed directed flow at high (26).
Ii The Saturation Momentum and Unintegrated Gluon Distribution Functions
The basic object required to calculate physical observables in the -factorized (28) picture employed by the gluon saturation models is the unintegrated Gluon Distribution Function (uGDF), . In principle, the uGDFs possess a Bjorken dependence determined by nonlinear evolution equations of the CGC theory (18); (19); (20) and their dependence is fixed by a characteristic saturation momentum, . In the McLerran-Venugopalan approach (17) the gluon distribution is suppressed below the saturation scale compared to the perturbative form . The parameterization of the KLN model as used in (7) is similar to the following Lorentzian form of .
The QCD coupling, , is regulated at lower scales by imposing a maximum value . The constant is a parameter set to reproduce as seen in experiment. The transverse coordinate dependence is implicit in the saturation momentum determined numerically for each nucleus.
In the KLN model, the saturation scale is assumed to be proportional to the participant density in a nuclear collision.
In Eq. LABEL:SatKLN, and is the collinear momentum fraction given by kinematics. The participant density is defined in terms of the usual Glauber thickness functions as in (7); (21). The integrated Gluon Distribution Function (iGDF) is calculated to be (in the perturbative regime),
In Eq. 3, the term accounts for the rapid growth of small Bjorken gluons while the factor of was introduced in KLN to account qualitatively for the rapid depletion of gluons as outside the small framework of the CGC model. The unintegrated gluon distribution shown in 1 is defined such that . This relationship holds in the leading logarithm DGLAP kinematic regime where . The KLN model has the freedom to set the normalization of independently of due to the factor in Eq. 3 in the leading log approximation. This variable is set so as to reproduce the experimentally implied value of . We will specify the numerical value of all the parameters after discussing the adjusted KLN type saturation model that we use in the paper. Our model differs from both the KLN and fKLN implementations of the CGC in this case, as we self consistently include the high depletion factor of in both and the determination of . Other models only include this factor in the uGDFs, leading to a relative depletion of particle production in the high regions in our model relative to previous ones. This feature will be essential to our description of asymmetric nuclear collisions later in this paper.
As can be seen in Eq. LABEL:SatKLN, the saturation momentum defined in the KLN model is not universal as it explicitly depends on the partner nucleus through its dependence on (23); (24); (26). The problem at hand is to modify the constitutive equations of the KLN model such that the saturation momentum depends only on the Glauber thickness function of the nucleus in question and still reproduces the phenomenologically successful participant dependence of the multiplicity that follows naturally from the conventional KLN model (21).
Just such an adjustment is detailed in (23) where they create the fKLN model. They show that the correct event averaging of the uGDF involved in a nuclear collision leads to the replacements and (using the notation found in (23)) where is the probability of finding at least one excited nucleon at the position in nucleus . The inelastic nucleon-nucleon cross section is set to mb. Explicitly writing out the changes,
Note that we no longer have the problems occurring from the limit that were found in (26) because the , or alternatively the , limit leads to a finite saturation momentum equal to that of a single nucleon.
Fig. 1 shows the squared saturation momentum as a function of the transverse coordinate at given value of Bjorken scaling variable, . One can easily see that as one nears the nuclear edge, the nuclear saturation momentum tends towards the nucleon saturation momentum from above. An added nice feature of the fKLN approach is that the nucleon saturation momentum is set uniquely by the low density limit of the nucleus saturation momentum and is not a free parameter as in (21).
Iii Multiplicity in Symmetric Collisions
We first concentrate on nucleon-nucleon collision and can calculate the multiplicity of charged particles in p-p collisions using the formula of Gribov-Levin-Ryskin (28).
The value of is set to reproduce data from the UA5 collaboration (29). One subtlety in the analysis is the role played by the transformation from the rapidity to the pseudorapidity . The conventional Jacobian is used as in (21) but there is the matter of setting the parameter found in the Jacobian. We use the experimentally determined average transverse momentum in the collision as a measure for (30) and find that the best fit is found when GeV. We use the experimentally found average (30) in the Jacobian for all collisions at GeV and assume that the ratio does not change as one goes to higher energies. The calculation for is shown alongside in the figured to show the effect of the Jacobian.
Fig. 2 shows the results of evaluating for different COM energies. One can see that the multiplicity is described well at GeV except at very high rapidities. This will be a recurring feature of the calculation as large rapidites are in the fragmentation region of the collision and are not well modeled by gluon saturation physics. The theory has been fit to the UA5 data (29) and all other calculations are predictions as there are no more parameters to tune. The predictions for the p-p collisions at LHC energies can also be seen in Fig. 2.
Once the parameters of the theory have been specified by reproducing p-p data, we move on to consider symmetric nuclear collisions; Au-Au collisions at the max RHIC energy of AGeV and Pb-Pb collisions at the expected LHC energy of AGeV. When we move over to a nuclear collision, we must treat particle production locally as a density in the plane transverse to the beam. The parameter values, in particular the normalization of the uGDF, is set to produce the correct number of particles in a p-p collision. In order to interpret this as a local transverse plane density, we note that the particle production in Eq. 6 is implicitly spread out over a transverse area the size of the nucleon-nucleon cross section, . We can further understand the smearing of particle production over an area of by noting (as shown in Eqs. LABEL:SatKLN and 4) that all nucleons present in an interaction region of size are considered part of a local density that determines the for particle production. We convert Eq. 6 into a local transverse space density,
The generalization of this to a full nuclear A+B collision is immediate.
We regard Eq. 8 as the defining expression for our extended local factorized KLN (elfKLN) model of the CGC. It includes the consistent local limit implemented by fKLN in addition to the self consistent high depletion factor of included in both the uGDFs and determination. This change makes it possible to set a universal normalization of the uGDFs such that the adjustments in going from to are twofold; we replace with and introduce a factor of in the nuclear distribution.
The left panel in Fig. 3 shows the predictions for and in Au-Au collisions at maximum RHIC energies along with data from the BRAHMS Collaboration (31) as a comparison. One can see that we achieve a very good description of the data over all centrality classes with the exception of the regions near the projectile and target rapidity where fragmentation effects dominate. The centrality cuts are imposed via a geometric Glauber analysis as in (21). We do not take into account the effects of the event by event fluctuation in geometry caused by considering the Monte Carlo Glauber analysis which has been shown to be important in explaining system size dependence of observables for smaller systems (especially eccentricity and elliptic flow with fluctuations) (32); (33); (34).
The right panel in Fig. 3 shows more clearly the system size dependence of particle productions at RHIC. It shows the charged particle production per participant pair as a function of for a number of different pseudorapidities. We get a good description of the particle production, as seen by the agreement between our calculations and data from BRAHMS (31) and PHOBOS (35). Note that, unlike previous attempts at reaching the limit (21), there is no problem in this calculation due to the consistent treatment of the local . The calculation extrapolates right to the charged multiplicity in a single p-p collision as one reduces . This agreement is non trivial and is a good check of our approach as all parameters have already been fit and self consistency demands that we get the correct limit.
The left panel of Fig. 4 shows the rapidity and pseudorapidity densities calculated from the integral of Eq. 8 in the case of a Pb-Pb collision at AGeV. The centrality cuts are performed using optical Glauber models as in (21). The predictions for LHC multiplicity are surprisingly lower than some previous predictions from saturation physics (21), even the ones using an fKLN type model (25). The cause of this difference can be tracked to the form for Eq. 8 imposed upon us by the self consistency we require between p-p and A-A collisions. Since we have a purely local description of particle production (which (21) does not include), we are forced to include the factor of to maintain consistency between Eqs. 6 and 8. Therefore, we are sensitive to two processes that increase total particle production as one increases COM energy; 1) the increase of gluon production per unit area caused by small evolution, and 2) the increase of the underlying interaction cross section used as an “averaging area” for coherent interaction (as can be seen by the presence of in the definition of the local density ). Previous calculations include 1) but not 2) and are thus typically larger by a factor of 70/42 1.66 than the current predictions.
The reason for this difference is schematically shown in Fig. 5. As one increases , one naturally lowers the typical value being probed in the interaction and thus increases the typical saturation scale. This leads to more particle production at higher energies as compared to lower energies. The coherence effect, however, is implicitly spread out over an area the size of . As schematically shown in Fig. 5, this increase in the coherence region limits the number of coherent production points one can “fit” into the interactions region, thus leading to a suppression in overall particle production.
Our elfKLN model consistently includes both effects mentioned. We systematically isolate the production per unit area by using Eq. 8 and thus can consistently explain both p-p and A-A data at RHIC energies, which previous approaches are unable to do without invoking an arbitrary enhancement for their p-p predictions (21); (25). The tendency is similar to the one found at RHIC and we note the consistent p-p system limit as we saw in the RHIC calculations. The right panel shows the particle production per participant pair at LHC energies as a function of .
Iv Multiplicity in Asymmetric Collisions
Eq. 8 can be directly generalized to the case of asymmetric collisions by inserting the Glauber thickness functions of the different collision partners used. For the particular case of deuteron-Au collisions at maximum RHIC energies, one needs to use a model of the nuclear density associated with the deuteron nucleus. We use the Hulthen distribution to describe the deuteron, using both and wave contributions as in (21). Once the distribution is determined, we calculate the Glauber thickness function and proceed as before.
Fig. 6 shows the comparison of our calculation with data from the PHOBOS Collaboration (38) for a number of different centrality cuts as well as for the minimum bias data. The first thing to notice is that there is very little agreement in the total number of particles produced at any centrality. The theoretical calculation consistently under predicts the data, so much as predicting a multiplicity for a centrality of that is below the multiplicity for a p-p collision. An interesting point to note, however, is that the shape of the overall data is very well reproduced just from the asymmetry of saturation momenta in the collision. The presence of the factor of in the determination of is essential to the reproduction of this asymmetry as it causes a precipitous drop in at high rapidities. Previous models (21) achieve asymmetric results by supplementing their far too slowly varying determinations with a leading logarithm evaluation for the integral in Eq. 8. This picks the points in phase space that possess the most asymmetric values of between the two collision partners.
We note that the determination of centrality cuts in our calculations was performed using an optical Glauber analysis as in (21). As previously shown in (21); (34) the deviations between a mean field treatment as in the current paper and a fully fluctuating Monte Carlo treatment as performed by the experiment is significant, specially when one is studying smaller systems. One way to determine the difference between the two approaches is to tabulate the average number of participants calculated in each, as shown in the table above. The experimental data in the table is taken from the PHOBOS Collaboration (38). One can see that there is an appreciable difference in the participant number between the two approaches, with the fluctuating approach isolating larger systems (i.e. smaller impact parameters) for each centrality cut. The table also explains the peculiarity found in Fig 6 for the centrality cut. The optical Glauber actually gives total participant number of less than 1 for that centrality cut, thus leading to particle production below the p-p level.
Given the different determinations of centrality in the approaches, a good way to compare the data to the calculation is to divide out by the overall number of participant pairs. This controls for the difference in the underlying impact parameter differences and poses the problem in terms of particle production per participant pair. The analysis is presented in Fig. 7 and it can be seen that the comparison is a lot more favorable now that differences due to centrality determination have been removed. There is good agreement between data and theory, specially in the mid to forward rapidity region. There is still some discrepancy as one moves backward in rapidity towards the nucleus. Intranuclear cascading increases the multiplicity of charged particle in the nucleus fragmentation region. Furthermore, we have neglected to take into account the geometric scaling region of the uGDF at low (39) which has been shown to be important in d-Au collisions at RHIC, at least off mid rapidity (40). Both these effects are larger for more central collisions, a fact that is borne out in Fig. 7 where the discrepancy between data and theory is reduced as one moves towards more peripheral collisions.
Having achieved a good description of d-Au data at RHIC we move on to a prediction for p-Pb collisions at LHC energies. Eq. 6 easily generalizes to an asymmetric p-Pb collisions by replacing one of the uGDFs with one that uses the saturation momentum of a nucleus. Fig. 8 shows predictions of a p-Pb collision for different centralities at GeV. The overall charged particle production is smaller by a factor of about 2 from p-Lb predictions made using the conventional KLN model (21). One must note, however, that this calculation makes centrality cuts using an optical Glauber model calculation and thus will run into similar problems to the d-Au calculations when being compared to data. The problems can be seen explicitly in the the centrality bin in Fig. 8, where the particle production is at a level lower than p-p at GeV (Fig. 2). The minimum bias prediction is most robust, as was the case with the d-Au predictions. We present this prediction as a good estimate of the expected multiplicity, one that will be better compared to data when presented as particle production per participant pair. In general, a complete description of asymmetric nuclear collisions may best be attempted using a model that can better accommodate fully fluctuating geometries, such as the Monte Carlo KLN (MCKLN) model detailed in (23).
The recent claims supporting the sQGP paradigm at RHIC are heavily dependent on the success of models using ideal hydrodynamics to model the evolution of the QGP stage of bulk evolution at RHIC (5); (6). This success is inherently linked to the choice of initial conditions assumed in the models, and an initial state derived from gluon saturation / CGC models does not work well in conjunction with ideal hydrodynamics to explain RHIC data (8). This is due to the larger initial spatial eccentricity produced in gluon saturation models which is naturally converted into a larger elliptic flow parameter, (22). A better understanding of the origin and robustness of this large eccentricity is needed. In order to do this we need to understand the and limits of the theory, which is equivalent to creating a model that works in describing not just central A-A collisions but also p-p and p-A collisions in one consistent model. This is especially important as recent work indicates that models that have more realistic nuclear edge characteristics lead to a spatial eccentricity that lies somewhere between the pure Glauber/BGK and KLN model eccentricities (23); (24).
We use the elfKLN model, our extension of the fKLN model in (23), consistently applied to p-p, p-A and d-A collisions along with the usual description of A-A collisions. We find that the model does an excellent job consistently describing particle production in symmetric nuclear collisions from p-p all the way up to A-A collisions. We fix our parameters to reproduce particle production in p-p collisions and predict the multiplicity in Au-Au collisions at RHIC, finding good agreement with the data. The model does very well in treating the nuclear edges consistently, as evidence by Fig. 3 where we show a consistency check between our Au-Au and p-p calculations. We also predict the multiplicity for Pb-Pb collisions at LHC energies, which are found to be surprisingly small. This is due to the increasing “coherence area” of the interaction, which we consistently build into our model.
We move on to calculate asymmetric nuclear collisions (d-A, p-A) in an effort to correctly model the edges of peripheral A-A collisions. We find that absolute agreement with the data is difficult to attain if one uses optical mean field methods to obtain the centrality cuts. However, good agreement is achieved with the d-Au RHIC data if one accounts for the variation in the underlying geometry by comparing particle production per participant pair.
In closing, we propose the elfKLN model as a robust model for bulk CGC physics that has the properties of being explicitly factorized and well defined in the diffuse nucleus limit. We do, however, recognize that it it might not be applicable in highly asymmetric (and thus smaller system size) nuclear collisions, where a fully fluctuating model such as the MCKLN of (23) might be more appropriate. Our prediction does not include effects such as running coupling and pre asymptotic terms in the evolution of the saturation momentum with energy (36); (37). These adjustments are left for future work.
Acknowledgements.Discussions with A. Dumitru, H.-J. Drescher, B. Cole, D. Molnar, A. Blaer, D. Winter and J. Jia are gratefully acknowledged. MG acknowledges support from the DFG, FIAS and ITP of the J.W. Goethe Uni. Frankfurt, and from GSI. This work is also supported in part by the United States Department of Energy under Grants No. DE-FG02-93ER40764.
- P. R. Sorensen, arXiv:nucl-ex/0309003.
- J. Adams et al. [STAR Collaboration], Phys. Rev. Lett. 92, 052302 (2004) [arXiv:nucl-ex/0306007].
- S. S. Adler et al. [PHENIX Collaboration], Phys. Rev. Lett. 91, 182301 (2003) [arXiv:nucl-ex/0305013].
- K. Adcox et al. [PHENIX Collab.], Nucl. Phys. A, 757 pp 184-283 (2005); B. B. Back et al. [PHOBOS Collab.], Nucl. Phys. A 757 pp 28-101 (2005); I. Arsene et al. [BRAHMS Collab.], Nucl. Phys. A 757 pp 1-27 (2005); J. Adams et al. [STAR Collab.], Nucl. Phys. A 757 pp 102-183 (2005);
- P. Huovinen, P. F. Kolb, U. W. Heinz, P. V. Ruuskanen and S. A. Voloshin, Phys. Lett. B 503, 58 (2001) [arXiv:hep-ph/0101136].
- U. W. Heinz, J. Phys. G 31, S717 (2005) [arXiv:nucl-th/0412094].
- T. Hirano and Y. Nara, Nucl. Phys. A 743, 305 (2004) [arXiv:nucl-th/0404039].
- T. Hirano, U. W. Heinz, D. Kharzeev, R. Lacey and Y. Nara, Phys. Lett. B 636, 299 (2006) [arXiv:nucl-th/0511046].
- D. Molnar and M. Gyulassy, Nucl. Phys. A 697, 495 (2002) [Erratum-ibid. A 703, 893 (2002)] [arXiv:nucl-th/0104073].
- J. M. Maldacena, Adv. Theor. Math. Phys. 2, 231 (1998) [Int. J. Theor. Phys. 38, 1113 (1999)] [arXiv:hep-th/9711200].
- G. Policastro, D. T. Son and A. O. Starinets, Phys. Rev. Lett. 87, 081601 (2001) [arXiv:hep-th/0104066].
- G. Policastro, D. T. Son and A. O. Starinets, JHEP 0209, 043 (2002) [arXiv:hep-th/0205052].
- G. Policastro, D. T. Son and A. O. Starinets, JHEP 0212, 054 (2002) [arXiv:hep-th/0210220].
- S. J. Brodsky, J. F. Gunion and J. H. Kuhn, Phys. Rev. Lett. 39, 1120 (1977).
- A. Adil and M. Gyulassy, Phys. Rev. C 72, 034907 (2005) [arXiv:nucl-th/0505004].
- X. N. Wang and M. Gyulassy, multiple jet production in p p, p A and A A Phys. Rev. D 44, 3501 (1991); Comput. Phys. Commun. 83, 307 (1994) [arXiv:nucl-th/9502021].
- L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 2233 (1994); 49, 3352 (1994); 50, 2225 (1994).
- E. Iancu, A. Leonidov and L. D. McLerran, Nucl. Phys. A 692, 583 (2001) [arXiv:hep-ph/0011241]. J. P. Blaizot and F. Gelis, Nucl. Phys. A 750, 148 (2005) [arXiv:hep-ph/0405305].
- I. Balitsky, Nucl. Phys. B 463, 99 (1996); Y. V. Kovchekov, Phys. Rev. D 60, 034008 (1999).
- J. Jalilian-Marian, A. Kovner, L. D. McLerran and H. Weigert, Phys. Rev. D 55, 5414 (1997); J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, Nucl. Phys. B 504, 415 (1997), Phys. Rev. D 59 014014 (1999), Phys. Rev. D 59 034007 (1999), [Erratum-ibid, D 59, 099903 (1999)]; J. Jalilian-Marian, A. Kovner and H. Weigert, Phys. Rev. D 59 014015 (1999); A. Kovner, J. G. Milhano and H. Weigert, Phys. Rev. D, 62, 114005 (2000); A. Kovner and J. G. Milhano, Phys. Rev. D, 61, 014012 (2000); H. Weigert, Nucl. Phys. A 703, 823 (2002); E. Iancu, K. Itakura and L. D. McLerran, Nucl. Phys. A 692, 582 (2001); E. Ferrero, E. Iancu, A. Leonidov and L. D. McLerran, Nucl. Phys. A 703, 489 (2002).
- D. Kharzeev, E. Levin and M. Nardi, Nucl. Phys. A 747, 609 (2005) [arXiv:hep-ph/0408050]. D. Kharzeev, E. Levin and M. Nardi, Nucl. Phys. A 730, 448 (2004) [Erratum-ibid. A 743, 329 (2004)] [arXiv:hep-ph/0212316]. D. Kharzeev, E. Levin and M. Nardi, Phys. Rev. C 71, 054903 (2005) [arXiv:hep-ph/0111315].
- A. Adil, H. J. Drescher, A. Dumitru, A. Hayashigaki and Y. Nara, Phys. Rev. C 74, 044905 (2006) [arXiv:nucl-th/0605012].
- H. J. Drescher and Y. Nara, Phys. Rev. C 75, 034905 (2007) [arXiv:nucl-th/0611017].
- T. Lappi and R. Venugopalan, Phys. Rev. C 74, 054905 (2006) [arXiv:nucl-th/0609021].
- H. J. Drescher and Y. Nara, arXiv:0707.0249 [nucl-th].
- A. Adil, M. Gyulassy and T. Hirano, Phys. Rev. D 73, 074006 (2006) [arXiv:nucl-th/0509064].
- B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. C 74, 021902 (2006).
- L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100, 1 (1983); Phys. Lett. B 100, 173 (1981); E Laenen and E. M. Levin, Ann. Rev. Nucl. Part. Sci. 44, 199 (1994).
- G. J. Alner et al. [UA5 Collaboration], Z. Phys. C 33, 1 (1986).
- J. Adams et al. [STAR Collaboration], Phys. Rev. C 70, 054907 (2004) [arXiv:nucl-ex/0407003].
- I. G. Bearden et al. [BRAHMS Collaboration], Phys. Rev. Lett. 88, 202301 (2002) [arXiv:nucl-ex/0112001].
- B. Alver et al. [PHOBOS Collaboration], J. Phys. G 34, S907 (2007) [arXiv:nucl-ex/0701049].
- P. Sorensen [STAR Collaboration], arXiv:nucl-ex/0612021.
- B. Alver et al. [PHOBOS Collaboration], Phys. Rev. Lett. 98, 242302 (2007) [arXiv:nucl-ex/0610037].
- B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. C 65, 061901 (2002) [arXiv:nucl-ex/0201005].
- D. Kharzeev, E. Levin and M. Nardi, arXiv:0707.0811 [hep-ph].
- J. L. Albacete, arXiv:0706.1251 [hep-ph].
- B. B. Back et al. [PHOBOS Collaboration], Phys. Rev. C 72, 031901 (2005) [arXiv:nucl-ex/0409021].
- E. Levin and K. Tuchin, Nucl. Phys. A 691, 779 (2001) [arXiv:hep-ph/0012167]. A. M. Stasto, K. J. Golec-Biernat and J. Kwiecinski, Phys. Rev. Lett. 86, 596 (2001) [arXiv:hep-ph/0007192]. E. Iancu, K. Itakura and L. McLerran, Nucl. Phys. A 708, 327 (2002) [arXiv:hep-ph/0203137]. A. Dumitru, A. Hayashigaki and J. Jalilian-Marian, Nucl. Phys. A 770, 57 (2006) [arXiv:hep-ph/0512129]. A. H. Mueller and D. N. Triantafyllopoulos, Nucl. Phys. B 640, 331 (2002) [arXiv:hep-ph/0205167].
- A. Dumitru, A. Hayashigaki and J. Jalilian-Marian, Nucl. Phys. A 765, 464 (2006) [arXiv:hep-ph/0506308].