QCD saturation at the LHC: comparisons of models to p+p and A+A data and predictions for p+Pb collisions

In a previous paper (arXiv:1011.1895), we showed that saturation models, constrained by e+p HERA data on inclusive and diffractive cross-sections, are in good agreement with p+p data at LHC in the soft sector. Particularly impressive was the agreement of saturation models with the multiplicity distribution as a function of $n_{\rm ch.}$. In this paper, we extend these studies further and consider the agreement of these models with data on bulk distributions in A+A collisions. We compare our results to data on central and forward particle production in d+Au collisions at RHIC and make predictions for inclusive distributions in p+Pb collisions at the LHC.

In a previous paper (arXiv:1011.1895), we showed that saturation models, constrained by e+p HERA data on inclusive and diffractive cross-sections, are in good agreement with p+p data at LHC in the soft sector. Particularly impressive was the agreement of saturation models with the multiplicity distribution as a function of n ch. . In this paper, we extend these studies further and consider the agreement of these models with data on bulk distributions in A+A collisions. We compare our results to data on central and forward particle production in d+Au collisions at RHIC and make predictions for inclusive distributions in p+Pb collisions at the LHC.

I. INTRODUCTION
In a previous paper [1], we discussed the computation of inclusive distributions in p+p collisions at RHIC and the LHC within the framework of k ⊥ factorization, wherein the unintegrated gluon distributions were determined from fits to small x HERA data on inclusive, diffractive and exclusive final states. The key ingredient in these fits is the dipole cross-section, which, to leading logarithmic accuracy, can be defined as whereŨ (b ⊥ ± r ⊥ 2 ) is a Wilson line in the fundamental representation representing the interaction between a quark and the color fields of the target. In the Color Glass Condensate (CGC) framework [2,3] of gluon saturation [4], the average · · · x is an average over these color fields; the energy dependence of the correlator as a function of x (or the rapidity Y = ln(1/x)) is given by the JIMWLK equation [5]. In the large N c limit, the equation for the energy evolution of this correlator is the Balitsky-Kovchegov (BK) equation [6].
We note however that neither JIMWLK nor BK is at present equipped to deal well with the impact parameter dependence of the dipole cross-section; the dipole cross-section in this formalism is taken in eq. (1) to be independent of the impact parameter. Another limitation of this framework is that the full next-to-leading logarithmic (NLL) expressions are not yet available; at the NLL level, only running coupling corrections to the leading log kernel have been considered in phenomenological applications. With these limitations in mind, we considered in ref. [1], saturation models of the dipole cross-section with the common criteria that their parameters be strongly constrained by fits to the HERA data. The saturation models considered 1 included the IP-Sat model [7], the b-CGC model [8][9][10] and the rcBK model [11]. All of these models provide good fits to the HERA data 2 .
In hadron-hadron collisions, one can derive at leading order the expression [15] This equation is a generalization of the well known k ⊥ factorization expression for inclusive gluon production [16] to include the impact parameter dependence of the unintegrated gluon distributions. Here C F = (N 2 c − 1)/2N c is the Casimir for the fundamental representation. Using a relation between quark and gluon dipole amplitudes strictly valid in the large N c limit, the unintegrated gluon distribution in either of the two protons can be expressed in terms of the corresponding dipole cross-section measured in DIS as [17] dφ(x, k ⊥ |s ⊥ ) 1 For more details, see the discussion in section II of ref. [1]. 2 The rcBK model has only been compared to inclusive HERA data [12]. We note however that this comparison is to the combined H1-ZEUS inclusive data [13], in contrast to the IP-Sat and b-CGC models, which were fit only to the older ZEUS data [14].
Thus the impact parameter dependent dipole cross-section determined from HERA data can be used to compute single inclusive gluon distributions in proton-proton collisions. Since this is a leading order computation, the overall normalization is not constrained and is determined from data as we shall describe below. For the integrated multiplicities, there is a logarithmic infrared divergence that can be regulated by introducing a mass term as discussed in our previous paper. We should mention that solutions of Yang-Mills equations that treat the infrared behavior properly give infrared finite distributions [18][19][20]. For a nice comparison of theoretical errors in various k ⊥ factorized approximations to the full classical Yang-Mills results, see ref. [21]. In ref. [1], we used eq. (2) combined with eqs. (3) and (1) to compute the rapidity and p ⊥ distributions for the models discussed for central rapidities in p+p collisions at RHIC energies all the way to the highest available LHC energies. A typical feature of the rapidity and p ⊥ distributions in these models was that the agreement with data improved with increasing energy; this should be the case because saturation effects are increasingly important at higher energies. We also observed that the saturation models, in particular the IP-Sat model, gave excellent fits to the multiplicity distribution P (n ch. ) as a function of n ch. for a wide range of energies. In this work, we will first briefly consider p+p collisions again before discussing A+A collisions and p/d+A collisions. In the latter case, we will present predictions for a future p+Pb run at the LHC.

II. RESULTS FOR P+P COLLISIONS
The probability distribution for producing n particles is where, in the Glasma flux tube framework [22], we can derive the negative binomial distribution [23] P NB n (n, k) = with the parameter k defined specifically to be with ζ=0.155 obtained from a fit to p+p multiplicity distribution. Here where, motivated by CGC computations on dilute-dense collisions [24], we choose Q S in the overlap area of the two hadrons to be Q S (s ⊥ , b ⊥ ) = min. {Q S (s ⊥ ), Q S (s ⊥ − b ⊥ )}. Also,n is the average multiplicity at a given impact parameter in the saturation model. Finally, as previously, we use dP dip. inel.
In fig. (1), we show results for the p+p multiplicity distribution plotted for √ s up to 7 TeV for |η| <0.5. This extends the comparison in ref. [1] to the CMS data at 7 TeV, with good agreement to n ch. = 60, nearly twice the range considered previously. The parameter ζ=0.155 is extracted from a fit to data at 0.9 TeV and used for all other energies 3 . It is the only free parameter in our fit; however, this quantity was recently computed non-perturbatively ab initio by solving Yang-Mills equations numerically [26] for two gluon correlations from gauge fields generated in the collision of two dense color sources. The results of the numerical computation vary depending on parameter choices in the range ζ ∼ 0.3-1.5-the lower end of this range is therefore a factor of two larger than the best fit value. Given the many uncertainties in the computation, the agreement is quite good, keeping in mind that nothing a priori prevents ζ from being orders of magnitude different. The good agreement of this framework with the LHC data over several decades, taken at face value, leads us to conjecture 4 that fluctuations in the number of produced gluons for a fixed distribution of hot spots gives a larger contribution to the multiplicity distribution than fluctuations in the distribution of hot spots themselves. The latter is of course only treated at the mean field level here in contrast to "pomeron loop" contributions [27] -suggesting perhaps that the latter are suppressed. Because Glasma flux tubes generate long range rapidity correlations, and can explain the distribution of high multiplicity events, our result provides further evidence corroborating computations in this framework [28,29] that suggest Glasma flux tubes generate the near side ridge seen in high multiplicity events by the CMS collaboration [30]. We now turn to transverse momentum distributions for charged hadrons and π 0 's, that in p+p collisions are computed from the expression where 5 D g→h (z, µ 2 ) is chosen to be 6.05 z −0.714 (1 − z) 2.92 . In ref. [1], we had considered the transverse momentum distributions only at mid-rapidity, where agreement of the IP-Sat and rcBK models was not particularly good. This agreement improved significantly at the higher LHC energies. The explanation provided there was that this better agreement is a consequence of smaller x values being probed at the higher energies. In fig. (2), we show the results of the rc BK and IP-Sat models for forward rapidities for RHIC p+p collision at 200 GeV. For both IP-Sat and rcBK we have used an overall normalization 6 extracted from energy dependence of the single inclusive multiplicity of the form: . This form absorbs the uncertainties in the inelastic cross-section and higher order effects (K-factors). For the IP-Sat (rcBK) model, one finds A = 0.23(6.15), b 0 = 5.77(5.14) and C = 0.32(0.76) using the mass term m = 0.4 GeV by fitting data points at η = 0 over the range of energy shown in fig. 3. For the rc-BK model the constant term A absorbs the prefactors 7 of Eq. 2 which includes unknown overlap area and other terms that cannot be separated from the " K factor ". The results are rather insensitive to the infrared cut-off m. However, one finds a ∼ 10% variation of the normalization when the 4 To confirm this conjecture, one would need to demonstrate that the results are valid for instance in the rcBK framework. However, because the impact parameter dependence of inclusive distributions in this model is rather simplistic and not constrained by data, this is difficult to confirm meaningfully at present. 5 The functional form quoted is for π + + π − ; we assume that the likelihood to fragment is identical for either charge and is equal to that for π 0 . For other charged hadrons, the functional form is assumed to be identical with the normalization constrained by the momentum sum rule. 6 The overall normalization for minimum bias p+p is slightly different from our previous paper [1] due to one less parameter (λ 0 ) in the large x extrapolation which is set to zero here. 7 Eq. 2 for min-bias p+p collision the rc-BK model normalization constant includes the constant prefactor ( 4α S πC F (2π) 5 ); here R A , R B corresponds to the radii of two protons and S A,B is the overlap area.
constants are extracted by a) fitting the full pseudo rapidity at RHIC energy, b) considering data points from different experiments. This variation, along with the numerical uncertainties, contributes to the gray bands shown in fig. 2. We see that the agreement at forward rapidities is significantly better than our previous comparison to the mid-rapidity distribution; this result provides a good benchmark for computing R pA at RHIC and in predictions of the same for p+Pb collisions at the LHC.
Transverse momentum distributions at forward rapidities in rcBK and IP-Sat models compared to STAR [34] and BRAHMS [35] data. The gray bands show the uncertainty in the determination of the normalization constant.

III. RESULTS FOR A+A COLLISIONS
For a large nucleus, in the IP-Sat model, we can approximate the dipole-nucleus cross section to be where σ dip (r ⊥ , x) p is obtained from integrating the dipole-proton cross section in eq. (1) over the impact parameter distribution in the proton. This form of the dipole-nucleus cross-section was shown previously [36] to give reasonable fits to the limited available fixed target e+A inclusive data. The initial conditions for rcBK evolution for a nucleus were similarly fixed by comparisons to the e+A data [37]. Substituting the expression for the dipole-nucleus cross-section in eqs. (2), and likewise the latter in (3), one can compute the nuclear multiplicity distributions. The infrared divergence in the multiplicity distribution is regulated in exactly the same was as was the case for the p+p multiplicity distribution, by replacing p ⊥ by m ⊥ = p 2 ⊥ + m 2 , with m varied between 0.2-0.4 GeV. Wherever we have considered fixed coupling, we have used α S =0.2; for the running coupling case, we run α S with the scale The results shown in fig. (3) are for 0-6% centrality, which corresponds to a median b ⊥ ≈ 12.2 GeV −1 ; we compute and N part (b ⊥ ) for this median value. We observe that a fairly good agreement with data is obtained for the infrared cut-off given by m = 0.4 GeV. The prescription for the running coupling gives a variation that corresponds to a 20% uncertainty at lower energies, which decreases significantly at higher energies.  [39][40][41][42] and for A+A from ref. [43,44] We next consider the centrality dependence of the multiplicity at RHIC ( √ s = 200 GeV) and LHC ( √ s = 2.76 TeV) in the IP-Sat model. While the agreement of the model with data shown in fig. (4) is reasonably good for the most central collisions, a systematic deviation is seen for lower centralities, and the model underpredicts the data. While within the range of the theoretical uncertainties outlined thus far, this systematic discrepancy leaves significant room for final state entropy production, which is expected to be more significant for more peripheral collisions. See for instance refs. [45,46] that estimate the amount of entropy production. As the right plot of fig. (4) shows, running coupling effects are less important for the most central collisions but introduce significant uncertainties relative to the fixed coupling results for more peripheral collisions.  We also note that a significantly better fit to the data at higher energies is obtained by including running coupling effects. We now employ eq. (4) to compute  [49,50] the multiplicity distribution in A+A collisions. While eq. (5) is computed identically to the p+p case, we need to determine the impact parameter distribution differently from the prescription used for the p+p case in eq. (7). The expression gives a better description of the impact parameter distribution in A+A collisions 9 . As in the computation of N part , σ NN ∼ 62 (41)  Interestingly, we find that the best fit is found for the value of ζ = 0.155 that also gives the best fit to the p + p data.
For the A+A collision data discussed here, we presented results for a) the energy dependence of central (0-6%) A+A collisions, b) the centrality dependence at two different energies, c) the rapidity dependence (at RHIC for 0-6% centrality and at LHC for 0-5% centrality), d) the multiplicity distributions (which are minimum bias). Unlike p+p collisions, for A+A collisions, we do not use an energy dependent normalization (because the nuclear size does not grow appreciably with energy unlike the proton) but a constant factor that is fit to the single inclusive distribution at one energy at η = 0 and used to fit these four data sets a) -d). For the IP-Sat model (the only model considered here), the value of this normalization (K-factor) is 0.07 for the mass term m=0.4 9 For A+A collisions, the published data [43] points are uncorrected requiring an additional parameter in contrast to the p+p case. The average multiplicity isn where the pre-factor Cm is the additional parameter specific to the multiplicity distribution and is tuned to provide a good fit to the uncorrected multiplicity distribution.

IV. RESULTS FOR p+A COLLISIONS
We computed the min-bias average multiplicity at mid-rapidity for the IP-Sat and rc BK models for p+A collisions 10 . The data is normalized to the PHOBOS 200 GeV d+Au data [51]. The energy dependence of the average multiplicity is shown in fig. (7), the band corresponds to the variation of m in the range of 0.2-0.4 GeV. Both the models give a comparable energy dependence, with the IP-Sat model giving a slightly higher multiplicity at the highest energies. The rapidity distributions for RHIC energies in the two models and predictions for LHC energies are shown in fig. (8). The models agree with the RHIC data with an accuracy of ≈ 10%, which is within the theoretical systematic uncertainty, which we shall discuss further shortly.  [51,52] Fig. (9) shows the transverse momentum distribution compared to BRAHMS [35] and STAR [34] 200 GeV data for h − and π 0 at forward rapidities. For the h − case, we have included 15% isospin correction in the normalization constant, to be discussed further below. Predictions for the p + A transverse momentum distributions for charged hadrons at η = 0 for LHC energies with this fixed normalization are shown in fig. (11). We compute the nuclear modification factor anticipated in p + A collisions at the LHC in the IP-Sat and rcBK models. R pA in both models, for √ s = 4.4 TeV/nucleon, approaches unity at p ⊥ ∼ 5 GeV. For √ s = 8.8 TeV/nucleon, the suppression persists to higher p ⊥ . The slope of R pA however appears quite different in the two models.
Regarding the overall normalization, for d+Au we have data available at only one energy. As noted, data and predictions include a) the energy dependence, b) the rapidity dependence at RHIC and LHC, and the the single inclusive pt distributions.
For the IP-Sat model the form of the normalization is A/(πb 2 max ) where, A = 0.25 with mass term m = 0.4 GeV. Unlike p + p here b max = 9.5 fm is a fixed number which is the maximum range of impact parameter to obtain minimum bias distribution and does not change with energy. For the rc-BK model a single normalization constant  A = 0.032 (with m = 0.4 GeV) absorbs 11 all the constant pre factors of eq.2. These normalization constants are obtained from a fit to the PHOBOS pseudo rapidity distribution; one obtains an ∼ 8% higher value of A when the fit is performed only to the data points at η = 0. Using the BRAHMS data for normalization also gives higher values for A which, along with other numerical uncertainties, contributes to the bands shown in fig. 9. The significant difference in A for IP-Sat and rcBK for d+Au collisions is because the the area of overlap and other terms in the pre-factor of kt-factorization are absorbed in A for rc-BK and cannot be separated from the "K factor". For the IP-Sat model, eq. (3) includes those factors and A is of the order of 1. This apparent difference in the two models doesn't affect any of our final results since same normalization is consistently used everywhere. Note also that in conversion from d+Au to p+A numbers we have used an additional factor of 1.6/2 in the normalization, which is standard in such conversions in the literature. In computing R pA , the result depends on N coll , which is sensitive to the proton inelastic cross-section. Since b max,proton grows with energy, one finds for the energies √ s = 4.4, 8.8 TeV that ratio b 2 max,proton (8.8TeV)/b 2 max,proton (4.4TeV) ∼ N coll (8.8TeV)/N coll (4.4TeV), a result consistent with expectations of N coll from Glauber approaches [53]; numbers quoted are in agreement with our ratio to 5%. 11 eq.2 for min-bias p+A collision the rc-BK model normalization constant includes the pre factor ( 4α S πC F (2π) 5 S p,A (πR 2 p )(πR 2 A ) ); here Rp, R A corresponds to the radii of the proton and nucleus and S p,A is the overlap area.

V. CAVEATS AND COMPARISONS
We would be remiss not to discuss the known sources of theoretical uncertainties in our framework; a brief discussion was also presented in ref. [1]. One that we alluded to previously is the k ⊥ factorization framework. Corrections to these come both from additional multiple scattering contributions in the dense projectile-dense target limit and from higher order contributions [54,55]; We estimate differences with the classical Yang-Mills multiple scattering effects to be no larger than other theoretical uncertainties and is accounted for an approximately 15% difference in normalization between the integrated multiplicity and the p ⊥ distributions. The higher order contributions are constrained by a normalization constant that's fit at a given energy. We have also shown in fig. (5) the effect of running coupling on multiplicity distributions. Another source of theoretical uncertainty arises from logarithmic sensitivity to the infrared cut-off seen in the single inclusive multiplicity (but not in the multiplicity distribution); the additional (classical Yang-Mills) scattering effect we alluded to also makes the distributions infrared finite. Even though the Yang-MIlls contribution is understood, it is at present cumbersome to include it in a "global" analysis.
In the IP-Sat model, as discussed previously in ref. [1], there is a sensitivity to the different data sets for the dipole cross-section that gave best fits to the then extant HERA data. We noted in ref. [1] that this sensitivity decreased with increasing energy. The sensitivity to the different fits to the gluon distribution is particularly severe at large p ⊥ where small differences are accentuated because of the steeply dropping cross-section. At high p ⊥ and large rapidities, the unintegrated distributions are sensitive to x ≥ x 0 (= 0.01). Guided by quark counting rule prescriptions, we choose 4 , as in eq. (20) of ref. [1], except with the parameter λ there set to zero. Further, as noted previously, the IP-Sat fits in ref. [9] were performed prior to the combined ZEUS-H1 data . One expects a revised fit to have an impact on the numbers extracted in ref. [9]. Finally, there is theoretical uncertainty arising from our imperfect knowledge of fragmentation functions.
Given the many uncertainties, one has to think of these predictions at this stage as having considerable elasticity of 20%-50% variability depending on the quantity studied. Only qualitative differences in measured distributions from theoretical projections (such as R pA > 1 at low p ⊥ at the LHC) can rule out models at present. Else, one has to check whether variations in parameters without adding new ones can accommodate the data. This is no different in spirit from global fits in perturbative QCD. Given that we have not tried to fine tune parameters, the agreement with the data presented here (in combination to the results in ref. [1]) is quite good considering that the dipole approach also gives a good description of small x HERA data.
We will now briefly compare our results with some recent related approaches. In ref. [56], results for the pseudorapidity distributions in p+p and A+A collisions at RHIC and the LHC are presented (based on previous work in ref. [57]) and predictions made for p+A collisions. The results shown are primarily for the b-CGC model which we also considered previously. While the b-CGC model does well for the pseudo-rapidity distributions, it does less well relative to the IP-Sat model for the multiplicity distribution P (n); as fig. (1) indicates, the latter works quite well. When it comes to nuclei, ref. [56] assumes that an additional gluon cascade mechanism to explain the faster energy dependence of the charged particle multiplicity. However, as pointed out in ref. [58], the difference in the energy dependence of the charged particle multiplicity in A+A collisions relative to p+p collisions is economically explained by the fact that the larger saturation scales in nuclei lead to more phase space for bremsstrahlung; while this may be conceptually similar to the cascade posited in ref. [56], no further modification of the formalism is necessary in the former approach. Indeed, we see in fig. (3), that a good agreement is obtained for both p+p and A+A data in the IP-Sat model. A similar quality of agreement 12 is seen in the rcBK model (combining results in refs. [1] and [59]). In ref. [60], results in the KLN model for the pseudo-rapidity distributions in p+p, d+A and A+A collisions at RHIC are presented, as are results for the LHC in p+p and A+A collisions, and predictions are made for the same in p+A collisions at the LHC. Similarly, results are presented for the P (n) in p+p collisions and predictions of the same made for p+A collisions. The KLN model gives good results for the pseudo-rapidity distributions at RHIC energies and in p+p and A+A distributions at LHC energies, which are similar to the models considered here. While the multiplicity distribution P (n) is motivated similarly to our discussion, the expression in ref. [60] does not depend on impact parameter; it will be interesting to see if a similar quality of agreement is obtained for higher n as in our fig. (1). The KLN model gives a better agreement to the centrality dependence of the charged particle multiplicity when compared to fig. (4); our result however leaves open the possibility that there may be significant entropy generation in the final state in more peripheral collisions. We also note that the KLN model does not at present allow for a global fit because it is not constrained by HERA data on e+p collisions [61].
Multiplicity distributions are sensitive only to p ⊥ < 1 GeV. For p ⊥ distributions, we first compare our results for R pA to those in ref. [62]. Only results for 8.8 TeV/nucleon are presented there, and for forward rapidities. Even so, we see that the suppression is considerably less in our framework. A significant fact is that in R pA , one has to account for a growth in the proton minimum bias cross-section with energy, which tends to suppress the denominator and lessen the suppression from what it would be otherwise. We found this effect to be quite significant and modify plots from those similar to ref. [62] to what we present here. We expect the same effect to impact the result (in fig. (5) of ref. [63] using a "hybrid" formalism [64] corrected by the inelastic contributions introduced in ref. [65]. The "hybrid" approach may be more suitable than our approach when distributions are sensitive to large x > 0.01 and high p ⊥ ; significant recent progress [55] in that approach will allow us to compute systematically corrections beyond the k ⊥ factorization framework. Finally, we observe that R pA in the IP-Sat model approaches unity more rapidly for √ s = 4.4 TeV/nucleon than for the leading twist shadowing evolved EPS08 curve [66,67] plotted in ref. [68], while the rcBK curve has a similar behavior to the EPS08 curve. This is not too surprising because shadowing in the IP-Sat model is from higher twist contributions that go away quickly at large Q 2 ; in contrast, there is a significant leading twist window in the rcBK model which may explain the qualitative similarity of the two results in this kinematic window.