Cronin effect vs. geometrical shadowing in d+Au collisions at RHIC

Multiple initial state parton interactions in p(d)+Au collisions are calculated in a Glauber-Eikonal formalism. The convolution of perturbative QCD parton-nucleon cross sections predicts naturally the competing pattern of low-pT suppression due geometrical shadowing, and a moderate-pT Cronin enhancement of hadron spectra. The formal equivalence to recent classical Yang-Mills calculations is demonstrated, but our approach is shown to be more general in the large x>0.01 domain because it automatically incorporates the finite kinematic constraints of both quark and gluon processes in the fragmentation regions, and accounts for the observed spectra in elementary pp-->\pi+X processes in the RHIC energy range, sqrt{s} = 20-200 GeV. The Glauber-Eikonal formalism can be used as a baseline to extract the magnitude of dynamical shadowing effects from the experimental data at differente centralities and pseudo-rapidities.

Multiple initial state parton interactions in p(d)+Au collisions are calculated in a Glauber-Eikonal formalism. The convolution of perturbative QCD parton-nucleon cross sections predicts naturally the competing pattern of low-pT suppression due geometrical shadowing, and a moderate-pT Cronin enhancement of hadron spectra. The formal equivalence to recent classical Yang-Mills calculations is demonstrated, but our approach is shown to be more general in the large x > 0.01 domain because it automatically incorporates the finite kinematic constraints of both quark and gluon processes in the fragmentation regions, and accounts for the observed spectra in elementary pp → πX processes in the RHIC energy range, √ s ∼ 20 − 200 GeV. The Glauber-Eikonal formalism can be used as a baseline to extract the magnitude of dynamical shadowing effects from the experimental data at differente centralities and pseudo-rapidities. It is well known that in proton (p), or deuteron (d), reactions involving heavy nuclei (A ∼ 200) at √ s < 40 AGeV, the moderate transverse momentum (p T ∼ 2 − 6 GeV) spectra are enhanced relative to linear extrapolation from p + p reactions. This Cronin effect [1,2,3] is generally attributed to multiple scattering of projectile partons propagating through the target nucleus. The data can be well accounted for phenomenologically by adding a random Gaussian transverse kick δk 2 T ∝ µ 2 L/λ to the projectile partons prior to hadronization [4,5,6,7]. Here λ is the parton mean free path in the nucleus, L ∝ A 1/3 is the average path length, and µ is a typical screening mass in ground state nuclei. These models naturally predict a slowly decreasing Cronin effect with increasing energy which has only recently been possible to test at √ s = 200 AGeV at the Relativistic Heavy Ion Collider (RHIC).
Interest in the Cronin effect has been revived, due to the developement of a new formulation of the physics based on the concept of gluon saturation and classical Yang-Mills field models [8,10]. In addition, a radical possibility was proposed in Ref. [9], that nonlinear gluon saturation may in fact strongly suppress moderate-p T spectra at RHIC. If such a deep gluon shadowing in this kinematic range were true, then an anti-Cronin suppression should have been observed in the large x = 2p T / √ s ≈ 0.01 − 0.1 and moderate 1 < Q = p T < 10 GeV scale range accessible at RHIC. Four experiments at RHIC [11,12,13,14] found independently that the nuclear modification factor, R dA (p T ) = 2dσ dA /Adσ pp , was consistent with a positive Cronin enhancement of hadrons with 2 < p T < 5 GeV. The magnitude of the enhancement is somewhat smaller than predicted for the π 0 , as reviewed in [15], but no evidence of strong shadowing was reported.
After the release of the RHIC data, the saturation model predictions of suppression were revised in Ref. [16,17,18]. In Ref. [16] Cronin enhancement with no high-p T suppression was shown to be a generic feature of saturation models. In Ref. [17] another version of saturation dynamics with Cronin enhancement coupled with the high-p T suppression of Ref. [9] was discussed. In Ref. [18], the Cronin enhancement at y = 0 was predicted to be progressively negated by non linear QCD evolution at smaller nuclear x, and therefore a gluon shadowing suppression is predicted at higher rapidities.
Different approaches to the calculation of the Cronin effect can be formulated in infinite momentum and target frames. In the traditional Glauber-Eikonal (GE) approach [19,20,21], sequential multiple partonic collisions in the target frame are computed. This leads to transverse diffusion and unitarity is naturally preserved. The low-p T spectra are suppressed by unitarity to compensate for the moderate p T Cronin enhancement. This is what we call "geometrical shadowing", as it is driven by the geometry of the collision. No high-p T shadowing is predicted in this approach.
In applications, the GE series has been directly evaluated thusfar only up to the three-scattering term and for √ s ≤ 40 GeV [3,19]. Numerically more convenient approximated GE models [4,5,6,7] have been proposed. They modify the pQCD rates through the inclusion of a nuclear broadened intrinsic k T , instead of evaluating the full GE series. Phenomenologically, it is well known [22,23] that intrinsic k T ∼ 1 GeV must be introduced to correct collinear factorized pQCD predictions to account for p + p data at moderate p T < 5 GeV. The approximated GE models simply extend that idea by adding a random kick δk 2 T to the intrinsic k 2 T . One drawback of such approaches is that a non-trivial p T or collision number dependence of the effective nuclear transport coefficient µ 2 /λ ∼ 0.05 GeV 2 /fm [5] must also be introduced to account for the actual Cronin data. While a logarithmic p T dependence of δk 2 T (p T ) is expected for partons undergoing multiple Yukawa screened interactions [21], the functional form of that p T dependence is usually adjusted to fit the Cronin data at one energy. A further drawback of such approximated GE models is that the unitarity constraints built into GE are ignored and hence the unitarity shadow and Cronin are treated as two seperate phenomena.
The more recent approaches [24,25,26,27] to the Cronin effect in the infinite momentum frame are based on the Mclerran-Venugopalan (MV) model of the nuclear wave funtion in classical Yang-Mills theory [28]. The general equivalence of GE and MV formulations for transverse diffusion was discussed in [30,31] in the context of gluon dominated small x ≪ 1 kinematics. In these approaches, the nucleus is approximated by a Weiszäcker-Williams gluon field with non-linearities approximated semi-analytically or computed numerically [27,29]. The non-linear gluon interactions lead to transverse diffusion and hence Cronin enhancement of nuclear partons prior to the scattering. The essential scale in this approach is a gluon saturation scale Q s = Q s (y, √ s, A), with y the rapidity of the produced gluon. One of the advantages of the MV approach is that unitarity is at least approximately enforced through the conservation of the number of virtual gluons in the transverse diffusion. Therefore these models predict a definite anti-Cronin suppression below some scale ∝ Q s . On the other hand, a disadvantage in present formulations of the MV model is that they ignore finite-energy kinematics of valence-and sea quarkinduced processes and the non-asymptotic large x > 0.01 features of gluon structure, where the classical approximation is unrealistic. A major disadvantage of MV models is that they cannot account for the elementary p + p transverse spectrum, that forms the denominator of the R pA nuclear modification factor. Neither can the models reproduce the absolute normalization of the spectra in pA collisions without extra phenomenological assumptions.
In this letter, we compute directly the GE series via numerical convolution of elementary parton-nucleon processes. An advantage over approximated GE models [4,5,6,7] is that our approach automatically conserves unitarity in the geometric optics sense of GE theory. In addition, we do not introduce extra phenomenolog-ical p T -dependent nuclear broadening of the intrinsic k T , since the GE series predicts the functional form of the Cronin enhancement based on the calculated nonasymptotic pQCD parton-nucleon cross sections. An advantage of our GE approach over MV model applications is that we autmatically include the finite kinematic largex feautures of both quark and gluon processes. Perhaps the most important advantage over the MV approaches is that our formulation is directly constrained to reproduce the absolute normalized spectra in p + p collisions. Therefore the GE approach presented below calculates consistently both p + p and p + A collisions at the finite energies accessible at RHIC.
Beside the geometrical quark and gluon shadowing, which is automatically included in GE models, at low enough x one expects genuine dynamical shadowing due to non-linear gluon interactions, as described in saturation models. Both kind of shadowing are present in the data, but it is not possible a priori to tell in which proprotion. The Glauber-Eikonal formalism can then be used as a baseline to extract the magnitude of dynamical shadowing effects from the experimental data.

PARTON-NUCLEUS COLLISIONS IN THE GLAUBER-EIKONAL MODEL
The GE expression for a parton nucleus scattering [20,21] is: where T A (b) is the target nucleus thickness function at impact parameter b. The differential and integrated parton-nucleon cross sections, dσ iN /d 2 k and σ iN (p 0 ) = dσ iN , are computed in pQCD as discussed below, see Eq. (5a)-(5b). The integrated parton-nucleon cross section depends on an infrared scale p 0 , which is determined by fitting p + p data. The exponential factor in Eq. (1) represents the probability that the parton suffered no semihard scatterings after the n-th one. In such a way, unitarity is explicitly implemented at the nuclear level, as discussed in Ref. [20]. The sum over n starts from n = 1 because we are interested in partons which are put on-shell by the interaction and later on hadronize.
The sum over n may be performed in Fourier space. The result reads: where, suppressing the dependence of the quantities in the r.h.s. on y, and It is possible to show [20] that, in the high-p T limit, Eq. (7) reduces to the usual single scattering approximation for parton-nucleus scattering: As p T →0 unitarity corrections switch on, suppressing the integrated parton yield [36], and inducing a random walk of the parton in p T space, thus redistributing the partons to higher p T compared to the single scattering approximation [20]. This is how the multiple scattering mechanism of Eq. (1) induces "geometrical" shadowing at low p T , and Cronin enhancement of the transverse spectrum at moderate p T , respectively. Note also that the accumulation of k T kicks in the multiple scattering process is computed in the model, not input as a Gaussian folding as in approximated GE models. The full expression for the p T spectrum, Eqs. (7) and (2a)-(2c), interpolates naturally between the geometrical shadowed low-p T and the Cronin enhanced moderate-p T regions. Note that σ iN (r T ) ∝ r 2 T as r T → 0 and σ iN (r T ) → σ iN as r T → ∞. This suggests the interpretation of σ iN (r T ) as a dipole-nucleon "hard" cross section. This dipole is of mathematical origin, and comes from the square of the scattering amplitude written in the Fourier variable r T , which represents the transverse size of the dipole. Then, we can interpret S iA , Eq. (2b), as the dipole-nucleus "hard" cross section, which clearly incorporates Glauber-Gribov multiple scatterings of the colour dipole. No other nuclear effects on PDF's are included beside multiple scatterings.
The interpretation of the interaction in terms of multiple scatterings of a dipole allows to relate this approach to other multiscattering formalisms such as Ref. [37,38] and the saturation computations of Ref. [18,25,26,31]. A first step in building such a dictionary was taken in Ref. [32], where it was shown that the dipole cross sections of Eqs. (2b)-(2c) are equivalent, in a suitable kinematic region, to the dipole cross section considered in the MV model of Ref. [25,26]. The main difference is in the input parton-nucleon cross section in p T -space. In our case, as we will discuss in the next sections, it is computed in pQCD, including full kinematics and interactions of the incoming parton with both quarks and gluons. Moreover, its energy and rapidity dependence are controlled by the DGLAP evolution of the parton distribution functions of the target nucleon. In the MV model, it is approximated using asymptotic kinematics for gluon targets only. Both models consider in S iA only inelastic dipole-nucleus scatterings. They neglect diffrective dipole-nucleus interactions, which however modify the p T -spectrum only at the lowest tarnsverse momenta [49].
Beside the geometrical quark and gluon shadowing, which is automatically included in GE models, at low enough x one expects genuine dynamical shadowing due to non-linear gluon interactions as described in the saturation models. However, it is difficult to disentangle these two sources of shadowing and suppression. The distinction between the two is however of fundamental interest as has already been emphasized in the context of e + p DIS HERA by A. Caldwell [44]. Most theoretical interest is not in the ubiquitous geometrical shadowing and unitarity corrections, but in the onset of genuine nonlinear QCD physics [45]. Moreover, saturation models cannot predict as yet the upper bound on x below which non linear effects set in. In order to help recognize possible novel nonlinear regimes it is essential to be able to calculate the baseline spectra isolating the unitarity and geometrical shadowing alone. The parton level GE model discussed below provides such a baseline.

INCLUSIVE MINIJET AND HADRON PRODUCTION IN pp COLLISIONS
Let's consider a pp ′ collision, where p and p ′ stand for a proton (p), a deuteron (d), or a nucleon (N ). In leading order pQCD, the inclusive cross section for production of a parton of flavour i = g, q,q (q = u, d, s, . . . ) with transverse momentum p T and rapidity y [39] may be written as a sum of contributions to the cross section coming from projectile (p) partons and from target (p ′ ) partons: Here we considered only elastic parton-parton subprocesses, which contribute to more than 98% of the cross section at midrapidity [34]. In Eq. (4), For the parton distributions we use the CTEQ5 parametrization at leading order [40]. The choice of the factorization scale Q p is discussed later. The cross sections dσ ij /dt of the ij→ij elastic partonic subprocesses can be found, e.g., in [23]. They are proportional to α s (µ 2 ), computed as in [39], at a scale µ = Q p The factor K in Eq. (4) is introduced in order to account for next-to-leading order (NLO) corrections [41], and is in general √ s and scale dependent [39,42].
Inclusive hadron production through independent fragmentation of the parton i into a hadron h, is computed as a convolution of the partonic cross section (4) with a fragmentation function D i→h (z, Q 2 h ): where q T is the transverse momentum of the hadron h, y h its rapidity, and z the light-cone fractional momentum of the hadron and of its parent parton i. For details, see Eqs. (8)-(11) of Ref. [39]. In this letter, we use LO Kniehl-Kramer-Pötter fragmentation functions [43] and set the fragmentation scale Q h = Q p . In the computation of the pp ′ cross section (6), we have two free parameters, p 0 and K, and a somewhat arbitrary choice of the factorization, renormalization and fragmentation scales. Our strategy is to compare two choices for those scales, namely Q p = Q h = m T /2 and Q p = Q h = m T , and then fit p 0 , K to hadron production data in pp collisions at the energy of interest. We analyze here π ± production at √ s = 27.4 GeV [2], and π 0 production at √ s = 200 GeV [35]. For the K-factor we perform a χ 2 fit to the high-p T tail of the data, following the procedure described in Ref. [39]. The fit of p 0 is performed by requiring that the computed spectrum does not exceed the experimental data at low p T 1 GeV. At √ s = 200 GeV, this fit is difficult because data exist for q T > 1.2 GeV only, so we used also data on charged hadron production [33]. The resulting data/theory ratio is plotted with thin lines in the bottom panel of Fig. 1, and the extracted parameters are listed in Table I. Note that the K-factor is strongly correlated to the choice of scale, while p 0 is more stable. Both of them depend on √ s.   I: Fitted values of p0 and the K-factor for π ± and π 0 production in pp collisions at √ s = 27.4 GeV and √ s = 200 GeV, respectively. We quoted the fit uncertainties only. The systematic uncertainty in the absolute normalization of experimental data (20% and 9.6%, respectively), which affect the determination of the K-factor are not included. The fit of the K-factor in the case of zero intrinsic momentum was made with the p0 determined using the k 2 T = 0.52 GeV 2 .
As in Ref. [39], we obtain a satisfactory description of data for q T 5 GeV over a broad range √ s, but the curvature of the hadron spectrum is overpredicted in the q T = 1 − 5 GeV range. As it is well known [22,23], this can be corrected by considering an intrisic transverse momentum k T for the colliding partons [42]. There exists many ways of implementing it phenomenologically, and we choose for simplicity a k T smearing of the cross section to approximate this effect. We introduce then unintegrated parton distributions where the width k 2 T of the Gaussian enters as a phenomenological parameter, and convolute over d 2 k 1T and d 2 k 2T in Eqs. (5a)-(5b).
We found that a fixed k 2 T = 0.52 GeV leads to a dramatic improvement in the computation of the trasverse spectra, which now agrees with data at the ±40% level. The quality of our pQCD computation including intrinsic k T is shown in Fig. 1, and the extracted K is reported in Table I. Without intrinsic k T it is not possible to fit the value of p 0 , due to the steepness of the data/theory ratio. The fit of the K-factor in this case was made with the p 0 determined using the intrinsic k T . The optimal choice of scale is found to be Q p = Q h = m T /2 at both energies, as the value of the K-factor is the closest to 1.
Let us comment briefly on the physical meaning of the infrared regulator p 0 . The divergence of the pQCD cross section for minijet production, indicates clearly a break-down of unitarity at low p T . A phenomenological way of restoring it, is to tame the divergence of dσ ij /dt  Table I. Bottom panels: the data to theory ratio for different choices of parameters. Solid lines are for Qp = Q h = m T /2, dashed lines Qp = Q h = m T . The pair of thin lines is computed with no k T smearing, k 2 T = 0 GeV 2 , and the pair of thick lines is computed with k 2 T = 0.52 GeV 2 . In the π 0 case [35], the dashed area shows the relative statistical and point-to-point systematic error added in quadrature, and does not include the systematic uncertainty of 9.6%, on the absolute normalization of the spectrum. In the π ± case [2], the shaded area includes statistical error only, without a systematic uncertainty of 20% on the absolute normalization of the spectrum.
by adding a small mass regulator, p 0 , to the exchanged transverse momentm p T . Seen in this light, p 0 represents the scale at which higher order parton processes enter into the game, beside the single scatterings considered in Eqs. (5a)-(5b). As the center of mass energy is increased, the partons probe the nucleon at smaller x, so that the density of target parton increases, and one should expect unitarity effects to arise at larger scales. The slight increase of the fitted p 0 = p 0 ( √ s) with energy is indeed a consistent check of this picture. For the same reason, and since quarks interact more weakly than gluons, one might expect also the unitarity corrections for quark-nucleon scattering to arise at p 0|quarks ≤ p 0|gluons . However, to check this relationship, we would need at least data on K ± production in a p T range of 1 − 6 GeV. Here we simply set p 0|quarks = p 0|gluons = p 0 .

FROM pp TO pA COLLISIONS
Having fixed the free parameters in pp collisions, we can proceed and compute the absolute transverse spectra in pA collisions. We assume the proton and the deuteron to interact as pointlike objects at an impact parameter b with the nucleus. Its nucleons, N , have isospin averaged parton distribution functions f i/N = Zf i/p +(A−Z)f i/n , with A and Z the atomic mass and atomic number. Furthermore, we assume that A-nucleus partons scatter only once on the proton or the deuteron, due to their small density. Then we may generalize Eq. (4) as follows, without introducing further free parameters: Hadron production is then computed analogously to Eq. (6). Note that, due to Eq. (3), at large p T , or as A→1 [assuming T A ( b)→δ( b)], the b-integrated pA cross section reproduces exactly the pp cross section discussed in the previous section. In this way, we can calculate consistently both the pp and pA transverse spectra in the same formalism.
The Cronin ratio, R BA , of the inclusive differential cross sections for proton scattering on two different targets, normalized to the respective atomic numbers A and B is given by First, we can test the GE formalism against low-energy data at √ s = 27.4 GeV [2] for the ratio of of midrapidity η = 0 pion spectra in proton-tungsten (pW ) and protonberillium (pBe) collsions:  [46], in order for a pA collision at fixed impact parameter to produce the same number of participant nucleons as a minimum bias one. The result is shown in Fig. 2. The two choices of scale, Q p = Q h = m T /2 and Q p = Q h = m T /2 -along with the respective fits of p 0 and K from Table I - tion for with b Au = 5.7 fm, to experimental data from the PHENIX collaboration [11]. The results obtained with the two choices of scales are similar, and only the computation with Q p = Q h = m T /2 is shown. At this energy the sensitivity of the result to the error in the fit of p 0 and to the scale choice is smaller than at Fermilab energy, thanks to the reduced steepness of the pp spectrum. The result is compatible with data on the whole p T range inside the experimental statistical and systematic errors. Despite this caveats, the GE model tends to slightly overestimate the data at p T 2 GeV. What we see in Fig. 3 is therefore a possible indication for a dynamical shadowing in addition to the basic Glauber geometrical shadowing. Its magnitude is consistent with the range of dynamical shadowing explored in [5,6,7] using a variety of shadowing functions [47].
In the right panel of Fig. 3, we plotted the corresponding Cronin effect at the parton level. The Cronin ratio peaks at fairly large p T ≃ 6 GeV, compatible with the expected z ≃ 0.6. The peak in our computation is positioned at significantly larger transverse momentum than found in the MV model of Ref. Ref. [26]. This difference may be due to their choice of parameters: they compute the Cronin ratio R BA for a nucleus A such that Q s/A = 2 − 3 GeV, and an arbitrary reference B such that Q s/B = 1 GeV. Also an infrared cutoff Λ = 200 MeV was employed. The same value of the peak (p T ≃ 3 GeV) was found in the MV model via numerical computation of Ref. [27], with a slightly different choice of parameters. None of the values of Q s/A , Q s/B and Λ where fixed by fitting absolute inclusive spectra in pp or pA collision. Therefore, the results of Ref. [26] should be understood as illustrative of the qualitative features of the MV model. As Fig. 3 clearly demonstrates, fragmentation strongly distorts the features of the parton level Cronin effect. Therefore, the transverse momentum scales illustrated in the qualitative saturation models as in Refs. [16,17,18,26], which up to now do not attempt to include the distortions of scales due to hadronization processes, should not yet be taken literally.
To understand better the possible emergence of dynamical shadowing at RHIC, we have two powerful handles. The first one is the rapidity dependence of the Cronin effect, which we will address in a separate publication (see also [6,18,26] for a discussion in the framework of approximated GE, MV and saturation models). The second handle is the centrality dependence of the Cronin effect: the more central the collision, the higher the density of target partons, the higher the shadowing effects induced by non-linear parton interactions. A very nice observable will be, in this repect, the ratio of p T spectra in central and peripheral collisions. This ratio has the additional experimental advantage that most of the systematic errors shown in Fig. 3 are expected to cancel out. This will provide a rather precise comparison to the GE model prediction which isolates geometric shadowing only. Our predictions is shown in Fig. 4 for two different choices of scale. The average impact parameter is b = 3.5 fm and b = 6.5 fm for central and peripheral collisions, corresponding to centrality classes 0-20% and 60-88%, respectively [48]. In the figure, the theoretical uncertainty due to the uncertainty in the fit of p 0 = 1.0 ± 0.1 GeV is shown as a dotted band. If dynamical shadowing is present, we would expect a larger deviation of the plotted curve from the data than what is observed in minimum bias collisions in Fig. 3.

CONCLUSIONS
We have studied the Cronin effect in pA and dA collisions in the context of Glauber-Eikonal models. These models incorporate parton multiple scatterings and unitarity in pQCD in a consistent way. Moreover, they include a detailed parton kinematics and reproduce, in the limit of A = 1, the hadron transverse spectra as computed in the pQCD parton model. The analysis of pp spectra allows to fix the free parameter of the model, and to compute the spectra in pA collisions without further assumptions.
A powerful feature of GE models is that they automatically include in the computation geometrical shadowing effects induced by unitarity and parton multiple scatterings. By isolating geometrical shadowing, one can use the GE model computations as a baseline in the search for genuine dynamical shadowing effects due to non-linear parton interactions.
We tested our computation of the Cronin effect in minimum bias collisions at √ s = 27.4 GeV. The same formalism applied to the recently measured d + Au data at √ s = 200 AGeV [11,12,13,14] describes well the Cronin effect at large p T , with a tendence to overestimate by ∼ 10 − 20% the effect for π 0 at η = 0 and p T 2 GeV. Our results are surprisingly similar to predictions based on phenomenological approximated GE models [4,5,6,7] in spite of the inclusion of geometrical shadowing in our GE approach. This provides further evidence for the possible existence of moderate shadowing in the x ∼ 0.01 − 0.1 range as explored in those references. However, radical gluon shadowing as predicted in [9] is not supported by the data. It remains to be seen if the most recent variations of saturation models can be fine tuned to account to the thusfar featureless R dAu ∼ 1 RHIC data. Future analysis of the centrality and pseudo-rapidity dependence of the Cronin effect at RHIC will provide a powerful tool to further constrain the magnitude of the dynamical shadowing effect.