Single Inclusive Hadron Production at RHIC and the LHC from the Color Glass Condensate

Using the unintegrated gluon distribution obtained from numerical simulations of the Balitsky-Kovchegov equation with running coupling, we obtain a very good description of RHIC data on single inclusive hadron production at forward rapidities in both p+p and d+Au collisions. No K-factors are needed for charged hadrons, whereas for pion production a rapidity independent K-factor of order 1/3 is needed. Extrapolating to LHC energies, we calculate nuclear modification factors for light hadrons in p+Pb collision, as well as the contribution of initial state effects to the suppression of the nuclear modification factor in Pb+Pb collisions.


Introduction
The suppression of particle production at forward rapidities in d+Au collisions compared to p+p collisions, experimentally observed at RHIC [1,2], constitutes one of the most compelling indications for the presence of non-linear QCD evolution effects in presently available data. The appropriate framework to study the nuclear wave function in this non-linear QCD saturation regime is the Color Glass Condensate (CGC), see e.g. the reviews [3,4] and references therein. The CGC is endowed with a set of non-linear pQCD evolution equations, the JIMWLK equations, which in the large-N c limit reduce to the Balitsky-Kovchegov (BK) equation [5,6]. The BK-JIMWLK equations can be interpreted as a renormalization group equation for the Bjorken-x evolution of the unintegrated gluon distribution, and more generally of n-point correlators averaged over the nuclear wave function, in which both linear radiative processes and non-linear recombination effects are included.
Indeed, the observed reduction of the forward hadron yield in d+Au collisions was predicted, albeit at a qualitative level, in [7,8], where it was directly related to the shadowing built in the wave function of the gold nucleus due to the enhanced role of non-linear effects in its evolution towards larger rapidities (smaller x). Later on, a better quantitative description of the d+Au forward hadron yields was achieved in the CGC calculations of [9,10,11,12]. These works relied on the use of models for the unintegrated gluon distribution of the gold nucleus inspired by approximate solutions of the BK equation. Relevant dynamical features in these models where either taken from analyses of lepton-proton scattering data or directly fitted to data, such as the anomalous dimension or the rapidity dependence of the saturation momentum Q s , the scale below which non-linear effects become important. More detailed analytical and phenomenological analyses of the corresponding nuclear modification factors were carried out in [13] and [14,15] respectively.
The reason why the BK-JIMWLK equations, the most robust theoretical tool available to describe the small-x dependence of the nuclear wave function, have not been directly used in phenomenological studies, is that they were originally derived at leading-logarithmic accuracy only. It was quickly understood that this was not good enough: in analytical [16,17,18,19] and numerical [20,21,22,23] studies of the original leading-order (LO) BK equation, the growth of the saturation scale was determined to be Q 2 s ∼ x −λ LO , with λ LO ≈ 4.8 Nc π α s , incompatible with the phenomenology of lepton-hadron scattering which demands Q 2 s ∼ x −0.2÷0.3 . Moreover there were hints that higher-order corrections would restore the compatibility of these values [24].
However, such insufficiency of the theory has been (at least partially) fixed through the recent calculation of the next-to-leading order (NLO) evolution equation. First, running-coupling corrections to the LO BK-JIMWLK kernel were derived in [25,26,27]. Then, the full NLO BK equation was obtained [28]. As demonstrated in [29], one can account for most of the higher-order contribution with running-coupling corrections only. Importantly, higher order corrections bring the evolution speed, λ, down to values compatible with experimental data, among other interesting dynamical effects, thus narrowing the gap between theory and data.
Indeed, first steps in promoting the BK equation with running-coupling corrections (referred to as rcBK henceforth) to an operational phenomenological tool have been taken in Refs. [30,31,32,33]. In [30], a good description of the rapidity and collision energy dependence of the hadron multiplicities in Au+Au collisions measured at RHIC was achieved. Ref. [31] demonstrated the ability of the rcBK equation to account for the small-x behavior of the total (F 2 ) and longitudinal (F L ) structure functions measured in e+p scattering experiments. Then it was shown in [32] that the proton scattering amplitude fitted to data in [31] allows a good simultaneous description of both the proton diffractive structure function measured at HERA and of the forward hadron yields measured in p+p collisions at RHIC. Finally, in [33] a good description of the few nuclear structure functions known at small x from e+A experiments was obtained.
Together, these works yield a consistent picture that present experiments can probe the nonlinear part of the hadronic and nuclear wave functions at small x, and that they can be successfully described by the CGC effective theory of QCD at high energies. In this work we provide a good description of the single inclusive hadron (charged hadron and neutral pions) yields measured in p+p and d+Au collisions at RHIC at forward rapidities (y h > 2), with unintegrated gluon distributions obtained from the rcBK equation. We also extrapolate our results to LHC energies and predict forward particle production in p+p, p+Pb collisions, that we present through nuclear modification factors for light hadrons. In the case of Pb+Pb collisions, we are able to give the contribution of initial state effects to the suppression of the nuclear modification factor.

Inclusive hadron spectra in d+Au collisions at RHIC
According to Ref. [34], the differential cross section for forward hadron production in protonnucleus collisions is given by where p t and y h are the transverse momentum and rapidity of the produced hadron, and f i/p and D h/i refer to the parton distribution function of the incoming proton and to the final-state hadron fragmentation function respectively. Here we will use the CTEQ6 NLO p.d.f's [35] and the DSS NLO fragmentation functions [36,37]. In writing Eq. (1) we have assumed that the factorization and fragmentation scales are both equal to the transverse momentum of the produced hadron. For light hadron production discussed here, the difference between the rapidity and pseudo-rapidity, η h , of the produced hadron can be neglected, yielding the following kinematics: describe the scattering of a hard valence quark (gluon) from the projectile on the saturated small-x glue of the target, either a nucleus or a proton. In order to avoid contamination from large(small)x effects in the target (projectile), we will restrict ourselves to the study of the forward region y h 2 both at RHIC and LHC energies, such that x 1 ≫ x 0 and x 2 ≪ x 0 , where x 0 is the xvalue where the small-x evolution starts (see below). Similar to previous approaches, we allow the possibility of a K-factor to absorb the effect of higher order corrections. For instance there is no α s -order term in Eq. (1), we shall only implement running-coupling corrections in the x 2 evolution ofÑ F (A) , but in principle they also affect the cross section [38].
The udg'sÑ F (A) are given by the two-dimensional Fourier transform of the imaginary part of the forward dipole-target scattering amplitude in the fundamental (F) or adjoint (A) representation, N F (A) , respectively: where r is the dipole size and Y is the evolution rapidity. In turn, the small-x dynamics of the dipole amplitudes is given by the rcBK equation: (3) For simplicity, we have omitted the subscripts F (A) in the r.h.s of Eq. (3). Using Balitsky's prescription [27], the kernel in Eq. (3) reads where r 2 = r − r 1 (throughout the paper we shall use notation v ≡ |v| for two-dimensional vectors).  Following [31], we regulate the running coupling in Eqs. (3) and (4) by freezing it to a constant value α f r s = 0.7 in the infrared. A detailed discussion about the different prescriptions proposed to define the running coupling kernel and of the numerical method to solve the rcBK equation can be found in [29]. The only piece of information left to fully complete all the ingredients in Eq. (1) are the initial conditions for the evolution of the dipole-nucleus(proton) amplitude. Similar to previous works, we take them from the McLerran-Venugopalan (MV) model [39]: where Q 2 s0 is the initial saturation scale (probed by quarks), and we take Λ = 0.241 GeV. Contrary to studies of e+p data, we have discarded initial conditions a la Golec-Biernat-Wüsthoff [40], since their Fourier transform result in an unphysical exponential fall-off of the ugd, and therefore of the hadron spectra as well, at large transverse momenta. Finally, in the large-N c limit which we have implicitly assumed in order to use the rcBK equation, the gluon dipole scattering amplitude can be expressed in terms of the quark amplitude as With this setup, we obtain a very good description of RHIC data. Fig. 1 shows the comparison of our results with data for the invariant yield of different hadron species in p+p and d+Au collisions at √ s N N = 200 GeV and rapidities y h = 2.2 and 3.2 for negative-charge hadrons (data by the BRAHMS collaboration [1]) and y h = 3.3, 3.8 and 4 for neutral pions (data by the STAR collaboration [2]). The only free parameters adjusted to the d+Au data are x 0 , the value of x which indicates the start of the small−x evolution, and Q s0 , the value of the saturation scale at x = x 0 . For the gold nucleus we obtain a quark saturation scale Q 2 s0 = 0.4 GeV 2 at x 0 = 0.02. Values of x 0 between 0.015 and 0.025 are allowed within error bands, they are used to generate the yellow uncertainty band in Fig. 1. A few comments are in order. First, the parameters Q s0 and x 0 are obtained from minimum-bias data, and therefore Q 2 s0 should be considered as an impact-parameter averaged value, the saturation scale at the center of the nucleus is bigger. We remind the reader that the corresponding gluon saturation scale is larger, Q 2,gluon s0 = 0.9 GeV 2 . Second, Q s0 and x 0 are compatible with other values extracted from e+A [33] or A+A [30] data. They can be compared with Q 2 s0 = 0.2 GeV 2 at x 0 = 0.007 obtained in the case of the proton (in [31], x 0 = 0.01 was obtained with Q 2 s0 = 0.2 GeV 2 ). Finally, no K-factor is needed in order to reproduce the charged hadron yields (i.e. K = 1), whereas a rapidity independent K-factor K = 1/3 is needed to describe the neutral pion data. Although the precise values of the K-factors do not have much meaning due to the 15% normalization uncertainties of the data, we do not have a good understanding of the strong hadron species dependence.

Nuclear modification factors in p+Pb and Pb+Pb collisions at the LHC
It is straigthforward to use formula Eq. (1) to calculate forward particle production in p+p and p+Pb at the LHC. We shall present our LHC results in terms of the nuclear modification factor where N coll is the number of binary proton-nucleon collisions in the p+Pb collision. In our predictions for p+Pb collisons at the LHC we use N coll = 3.6, which is half the number of collisions determined in minimum bias d+Au collisions at RHIC [1]. Thus, in order to compare our results with experimental data one should renormalize our curves in Figs (2) and (3)  number of collisions determined experimentally at the LHC. Note that computing the ratio Eq. (7) removes the sensitivity to the K factors. Our R pP b calculations are displayed for two possible LHC energies ( √ s N N = 8.8 and 6.2 TeV) in Fig. 2 for pion production and Fig. 3 for charged hadron production. In both cases, one observes the expected trends that R pP b decreases with increasing y h , and increases with increasing p t . Our results for R pP b indicate that a significant suppression ∼ 1/2 should be expected already at not too forward rapidities. It is highly likely that the CGC dynamics studied here would also lead to a similar suppression at mid-rapidity (see, e.g. [15]). However, to make a clear quantitative prediction for mid-rapidity one should ensure a proper treatment of high-x effects in the target, which is beyond the scope of this paper. We also compare the y = 2 and 4 curves with predictions obtained with the k t −factorization formalism, in order to check the validity of that approach, and especially to test up to what value of y it can be used. The k t −factorization formula (see Eq. (8) below) is valid when the dominant contributions to the cross section come from small values of x, for both the projectile (x 1 ≪ 1) and the target (x 2 ≪ 1). For instance, it only includes gluonic degrees of freedom. This approach is clearly insufficient at very forward rapidities or large p t , where valence quarks of the projectile are important (x 1 1). However, as can be seen in Figs. 2 and 3, both formalisms give comparable results, as the lines from k t -factorization overlap with the uncertainty bands spanned by the results from the hybrid formalism. This seems to identify a kinematical window where both approximations (Eq. (1) and Eq. (8) below supplemented with parton fragmentation) are valid. To some extent it is not surprising that both formalisms yield comparable nuclear modification factors, since the suppression is ultimately rooted in the udg's themselves.
The advantage of the k t -factorization formalism is that it can be used at mid rapidities, when x 1 also becomes very small, invalidating formula Eq. (1). This is true at the LHC where at mid-rapidity x 1 ∼ x 2 ≪ 1, but not at RHIC where both values of x are generally too large at y = 0. Indeed one should keep in mind that p t e y h / √ s N N and p t e −y h / √ s N N are only lower values for x 1 and x 2 respectively, but that through fragmentation larger values of x actually contribute more.
The k t -factorization formula to describe inclusive gluon production reads [41] dN AB→gX dη d where B t is the impact parameter of the collision. Gluon fragmentation is not explicitely written down (therefore x A = p t e y / √ s N N and x B = p t e −y / √ s N N ), but it should also be accounted for.
The unintegrated gluon distribution in Eq. (8), ϕ, is actually simply related to the one in Eq. (1): A detailed discussion about the definition and physical interpretation of the different udg's discussed here can be found in [42]. We had not specified the impact parameter dependence of the ugd's before because it is not needed in a p+p or p+A collision. Indeed in these cases one can write . The b-integrated proton udg does not appear in formula Eq. (1), it is rather the standard p.d.f.s that describe the (dilute) proton content in this formalism, while in the k tfactorization case the b dependence of ϕ p (b) can be safely neglected if we are only looking at ratios such as Eq. (7). As for the collision impact parameter B t dependence of the target ugd ϕ B , we have been dealing with minimum-bias data, therefore as mentioned before, the ugd's obtained with Q 2 s0 = 0.4 GeV 2 for the gold nucleus (and 0.2 GeV 2 for the proton) should be thought of B t averaged udg's, but in principle we could also look at different centrality bins. Note that formula Eq. (8) has only been proven to be valid for p+p and p+A collisons (or dilute-dense scattering) [43,44,41,45,46], and there are several hints that it does not hold for A+A collisions (or densedense scattering) [47], even if there were no final-state effects. However, numerical results seem to indicate that k t −factorization breaking may not be too important in practice [30,48,49], but this could be process dependent. With the above remarks in mind, we shall use Eq. (8) to calculate the initial state effects on particle production at mid-rapidity in Pb+Pb collisions at the LHC.
We deal with the impact parameter in the following way: we assume that it factorizes as This is acceptable for minimum bias results, and in this case the functionsφ A andφ p are simply related to the averaged ϕ A and ϕ p , used for minimun bias p+A and p+p collisions: N pA collφ A /φ p = ϕ A /ϕ p . We shall again use N pA coll = 3.6 in our LHC calculations. The nuclear modification factor R P bP b we obtain is shown in Fig. 4 at the gluon level, it corresponds to gluon production immediately after the collision, i.e at proper times τ = 0 + . Obviously, in order to compare our results with data, one should convolute them with final state effects due to interactions of the produced gluons with the Quark Gluon Plasma, and with hadronization effects as well. Although this is beyond the scope of this work, our results indicate that a sizable part of suppression expected for single hadron yields at the LHC [50] may be due to purely initial state effects. Finally, note that the value of the saturation scale probed by gluons at x 0 = 0.02 is Q 2 s0 = 0.9 GeV 2 , this is used in Fig. 4 for minimum-bias predictions, rather a band is generated using Q 2 s0 = 0.8 and 1 GeV 2 . To study different centrality bins in Pb+Pb collisions, first one would have to improve our approximation for the b dependence of the ugd's.

Conclusions
In this work we have presented a good description of the hadron yields measured at forward rapidities in p+p and d+Au collisions at RHIC using the hybrid formalism proposed in [34] to describe high-energy dilute-dense scattering. The main new ingredient in our calculation is the use of the BK equation including running-coupling corrections to describe the Bjorken-x dependence of the nuclear (proton) wave functions. With the two free parameters in our work, x 0 and Q s0 fitted to RHIC data, we extrapolate to LHC energies without further adjustments and predict the suppression of the different forward hadron yields in p+Pb collisions with respect to p+p collisions. Using a different formalism, k t -factorization, we estimate the contribution of initial state effects to the nuclear modification factor in Pb+Pb collisions at the level of gluon production at the LHC.
While our results offer an additional indication for the presence of CGC effects in RHIC data, the interest now focuses mostly in calibrating the expectations for the Heavy Ion program at the LHC. Our predictions for the LHC rely on the most up-to-date tools within the CGC formalism (Pomeron-loop corrections have been looked at [51], but are only relevant at asymptotically large energies when the running coupling is included [52,53]) and will be useful to confirm the tentative conclusions reached at RHIC and to distinguish between alternative physical scenarios, like those proposed in [54,55] (a complete set of predictions stemming from different formalisms can be found in [50]), where the suppression of forward yields at RHIC is due to the non-eikonal propagation of the leading parton, resulting in energy loss in the forward region. Finally, the proper characterization of gluon production in the early stages (before thermalization) of Pb+Pb collisions at the LHC would serve as crucial input for studies of the medium produced in such collisions, such as hydrodynamic simulations or studies of jet quenching. In the latter case, the suppression predicted here due to initial state effects would add to the final state effects due to the presence of a Quark Gluon Plasma.