Forward $J/\psi$ production in proton-nucleus collisions at high energy

Inclusive production of $J/\psi$ mesons, especially at forward rapidities, is an important probe of small-x gluons in protons and nuclei. In this paper we re-evaluate the production cross sections in the Color Glass Condensate framework, where the process is described by a large x gluon from the probe splitting into a quark pair and eikonally interacting with the target proton or nucleus. Using a standard collinear gluon distribution for the probe and an up to date dipole cross section fitted to HERA data to describe the target we achieve a rather good description of the cross section in proton-proton collisions, although with a rather large normalization uncertainty. More importantly, we show that generalizing the dipole cross section to nuclei in the Glauber approach results in a nuclear suppression of $J/\psi$ production that is much closer to the experimental data than claimed in previous literature.


I. INTRODUCTION
Heavy quarks occupy an important place in QCD phenomenology. The mass of the charm quark is large enough that one should be able to describe its interactions using weak coupling. On the other hand, it is light enough to be a dynamic probe that is sensitive to momentum scales that characterize high energy nuclear collisions, such as the deconfinement temperature or the gluon saturation scale. Because of the clean experimental signatures, vector meson production is a particularly important probe of heavy quark interactions. The "melting" of the J=ψ meson has traditionally been considered one of the most important signals of the deconfinement transition in ultrarelativistic heavy ion collisions. Studying J=ψ production in proton-proton and proton-nucleus collisions at high energies can also provide valuable insight into the physics of gluon saturation and strong color fields at small x, as well as into the energy loss of a high energy probe in nuclear matter.
We will in this paper use the color glass condensate (CGC) framework together with a simple color evaporation model (CEM) to calculate cross sections for J=ψ and ϒ production at forward rapidities in proton-proton and proton-nucleus collisions at LHC energies. At forward rapidity the probe proton can be treated as a "dilute" collection of partons, and the CGC calculation proceeds according to a rather straightforward picture. The probe proton can be described by conventional collinear parton distributions in the so-called "hybrid" formalism, or by a corresponding "k T -factorized" approximation (which in fact strictly speaking violates k T factorization [1] since it involves several unintegrated gluon distributions). The physical picture in both cases is the same: a gluon from the probe can split into a quark-antiquark pair. As this system propagates through the target, it can pick up an eikonal phase described by an SU(3) color matrix, the "Wilson line" from the target. This leads to an expression for the cc pair production cross section (in either singlet or octet color states) in terms of Wilson line correlators describing the target. The same Wilson line correlators, the simplest one of them being the "dipole cross section," appear in calculations of e.g. total deep inelastic scattering (DIS) cross sections [2][3][4], single [5][6][7][8] and double inclusive [9][10][11][12] particle production in proton-proton and proton-nucleus collisions, diffractive DIS [13,14] and initial state for the hydrodynamical modeling of a heavy ion collision [15][16][17], making the framework very broadly applicable to many different processes.
Several LHC experiments have performed detailed measurements of inclusive J=ψ production in a broad range in rapidity, particularly in proton-proton [18][19][20][21][22] but also in proton-nucleus [23,24] collisions. The results on the "cold nuclear matter suppression," i.e. the suppression of the cross section in proton-nucleus relative to protonproton collisions, have then been compared to several theoretical predictions. In particular, the nuclear suppression predicted by the CGC calculation in [25] has been much stronger than observed in the data. The main purpose of this paper is to explore to what extent this disagreement stems from the framework itself, and how much is due to approximations concerning nuclear geometry that are necessary for relating the CGC description of the proton to that of the nucleus. Indeed, in the case of single inclusive light hadron production [8], the disagreement with early CGC nuclear suppression calculations [26] and LHC data turned out to be mostly due to the nuclear geometry effects. We will in this paper calculate cross sections for inclusive (prompt) forward J=ψ production in proton-proton and proton-nucleus collisions using the collinear hybrid framework also used in Ref. [25].
We use here the dipole cross sections obtained in Ref. [8] by performing a fit to HERA data using the running coupling Balitsky-Kovchegov (BK) equation. These have several advantages over the ones used in previous works. First, the dipole cross section in a nucleus is related to the one in a proton taking into account the fluctuating positions of the nucleons in the nucleus via a Glauber model, without any additional parameters besides the standard Woods-Saxon nuclear density. This leads, unlike many previous CGC calculations, to nuclear modification ratios that explicitly approach one in the high transverse momentum limit, at all collision energies. Physically this should be expected from the validity of collinear factorized perturbative QCD in this regime. Second, the dipole cross section from Ref. [8] corresponds to an explicitly positive unintegrated gluon distribution when Fourier transformed to momentum space. Third, our calculation, like that of [8], takes consistently into account also the proton transverse area from the fit to HERA data. This completely fixes the normalization of the leading order (LO) calculation at the quark level without any additional parameters to describe the proton size. 1 The use also of the proton transverse area from the DIS fits and the way of extending the dipole cross section from a proton to a nucleus are quantitatively the most important differences between our work and that of [25].
This paper is organized as follows. We will first briefly describe the formalism used in the calculation in Sec. II and the used dipole correlators in Sec. III. We will then present our results for cross sections in proton-proton and protonnucleus collisions in Sec. IV and finish with a discussion and outlook for future work in this direction in Sec. V.

II. FORMALISM
We will here use the simple CEM, where a fixed fraction of all cc pairs produced below the D-meson threshold are assumed to become J=ψ mesons. The quark pair can be produced in either the singlet or octet representation, the additional color is assumed to be "evaporated" away in the form of a soft gluon. In the CEM model, the differential cross section for J=ψ production therefore reads where P ⊥ and Y are the transverse momentum and rapidity of the produced J=ψ, dσ cc d 2 P ⊥ dYdM 2 is the cross section for the production of a cc pair with transverse momentum P ⊥ , invariant mass M and rapidity Y, m c is the charm quark mass and m D ¼ 1.864 GeV is the D-meson mass. Here F J=ψ is a nonperturbative quantity representing the probability that a cc pair has to form a J=ψ. In the literature [27] this constant is taken to have different values depending on the c-quark mass, the parton distribution and the factorization scale. We will here use a constant value F J=ψ ¼ 0.06. This is somewhat larger than the typical values in perturbative NLO calculations, but it should be kept in mind that our calculation is at the LO level and is therefore likely to underestimate the normalization of the cc cross section. The same model can be used to compute the cross section for ϒ production by simply replacing m c with the bottom quark mass, M D with the B meson mass M B ¼ 5.280 GeV and F J=ψ with F ϒ , which we take equal to 0.04 here. One should note that since the b-quark mass is already significantly higher than the saturation scale, bottom quark production is not in the natural regime of validity of the CGC framework. We will in this paper however calculate also ϒ cross sections for comparison.
The formalism for analyzing gluon and quark pair production in the dilute-dense limit of the CGC formalism was worked out in great detail already some time ago [28,29] (see also [30]) and applied in many calculations such as [25,[31][32][33]. It can be obtained from the observation that an incoming gluon from the projectile can split into a quark-antiquark pair either before or after the interaction with the target. These partons are then assumed to eikonally interact with the target, picking up a Wilson line factor in either the adjoint or the fundamental representation, depending on the particle. We will here use the collinear approximation, where the incoming gluon is assumed to have zero transverse momentum. Calculated in this way the cross section for the production of a cc pair with transverse momenta p ⊥ and q ⊥ and rapidities y p and y q reads, in the large-N c limit [25], The longitudinal momentum fractions of the gluons coming from the projectile and the target, x 1 and x 2 , are given by The "hard matrix element" Ξ coll is the sum of three terms, corresponding to a quark pair scattering off the target (qq; qq), the gluon scattering before splitting into a quark pair (g; g) and an interference between the two. The explicit expressions are ð5Þ The propagation of a quark-antiquark pair through the color field of the target is described by the function where b ⊥ is the impact parameter. All the information about the target is contained in the function S Y ðk ⊥ Þ, which is the fundamental representation dipole correlator in the color field of the target: where Uðx ⊥ Þ is a fundamental representation Wilson line in the target color field. Our calculation is a leading order one, so for consistency we should use a leading order collinear gluon distribution to describe the probe. In this work we use the MSTW 2008 [34] LO parametrization for G p ðx 1 ; Q 2 Þ.

III. DIPOLE CORRELATOR
For the dipole correlator S Y ðk ⊥ Þ we use the MV e parametrization from Ref. [8]. This is obtained by fitting the combined inclusive HERA DIS cross section data [35] in the region Q 2 < 50 GeV 2 and x < 0.01 using a dipole cross section, whose x-dependence is obtained with the running coupling BK equation [36][37][38]. The dipole amplitude at the initial rapidity x 0 ¼ 0.01 is given by the parametrization The impact parameter profile in the proton is assumed to factorize from the dimensionless dipole amplitude, and is parametrized by a constant so that the distribution needed for quark pair production is given by Using the expression for the running coupling constant, the fit gives a very good description χ 2 =degrees of freedom ¼ 1.15 of the data with the parameters Q 2 s0 ¼ 0.060 GeV 2 , C 2 ¼ 7.2, e c ¼ 18.9 and σ 0 =2 ¼ 16.36 mb. One advantage of the parametrization [8] over the AAMQS [2] one used in [25] is that the Fourier transform of S Y¼ln 1 x 0 ðr ⊥ Þ is positive definite even at the initial rapidity, enabling a more natural physical interpretation as an unintegrated gluon distribution. We also used the McLerran-Venugopalan (MV) and MV γ parametrizations from Ref. [8] for the initial condition of the BK equation and found that for LHC energies (where the dipole correlator is already evolved over a significant rapidity interval from the initial condition) the effect of such a change is small compared with the other uncertainties involved in the calculation.
For the nucleus we use, again following [8], a BK-evolved dipole correlator with the initial condition where with d ¼ 0.54 fm and R A ¼ ð1.12A 1=3 − 0.86A −1=3 Þ fm, is the standard Woods-Saxon transverse thickness function of the nucleus, normalized to unity. To calculate cross sections with nuclear targets we explicitly integrate over the impact parameter, treating the dilute edge of the nucleus (where the saturation scale of the nucleus falls below the proton saturation scale) as in Ref. [8]. In practice this amounts to using the proton-proton result scaled such that the nuclear modification factor R pA , defined in Eq. (15), is unity in the dilute edge of the nucleus. This "Glauber model" form can be derived by assuming that at the initial rapidity, a high energy probe scatters independently off the nucleons in the nucleus. The average over the fluctuating positions of the nucleons can then be performed analytically. At small r the integral of Eq. (13) over b explicitly approaches A times the dipole-proton cross section (9). Thus the form (13) explicitly leads to nuclear modification ratios R pA [see definition Eq. (15)] that, at all collision energies, approach unity at small distances, i.e. large momenta. We emphasize that besides the standard nuclear density T A ðbÞ there are no additional parameters in passing from the proton, Eq. (9), to the nucleus in Eq. (13). Also the transverse area of the proton σ 0 =2 is consistently taken to be the same in proton DIS and in hadronic collisions, i.e. the identification (10) is also used in the pp → J=ψ þ X cross section. In our calculations the nuclear modification only depends on parameters fit to proton data and on the standard nuclear thickness function. In particular, no independent nuclear saturation scale with an unknown value is needed, in contrast to Ref. [25]. The practical result of this is that the effective typical nuclear saturation scale from Eq. (13) is smaller than the A 1=3 scaling from the proton assumed in Ref. [25]. On the other hand, we treat separately the area occupied by the small-x gluons in the proton, σ 0 =2, and the total inelastic cross section, σ inel . Thus here the J=ψ cross section in proton-proton collisions is proportional to σ 0 =2, which is smaller than the proton area assumed in Ref. [25]. The net effect of the smaller nuclear Q s and the smaller yield in proton-proton collisions is that one obtains a nuclear suppression ratio R pA that is closer to one around P ⊥ ∼ Q s than in Ref. [25], and simultaneously approaches unity at large P ⊥ .

IV. RESULTS
In this section we will present the results of our calculation and show a comparison with LHC data from the ALICE and LHCb experiments in pp and pA collisions at center of mass energies of 7 and 5 TeV respectively. It should be noted that both for pp and pA collisions, LHCb gives separately the cross section for prompt J=ψ production and the contribution from b decays (which is not included in our calculation) while ALICE data is only for inclusive J=ψ production. Therefore, for a consistent comparison, we subtract the b decay fraction measured by LHCb from the ALICE data to extract the prompt contribution.
In the uncertainty band we include the variation of the charm quark mass between 1.2 and 1.5 GeV (in the case of J=ψ production) or of the bottom quark mass between 4.5 and 4.8 GeV (in the case of ϒ production) and of the factorization scale between 2M ⊥ and M ⊥ =2 where M ⊥ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and M is the qq pair's invariant mass. When we show quantities integrated over the transverse momentum of the produced meson, we integrate up to P ⊥ ¼ 15 GeV. Our results do not depend significantly on this choice since the cross section is quickly decreasing at large P ⊥ .

A. Proton-proton collisions
We show our results for the P ⊥ -integrated cross section dσ J=ψ dY at a center of mass energy ffiffi ffi s p ¼ 7 TeV in Fig. 1, compared with data from ALICE [22] and LHCb [19]. We see that the shape of the data is quite well described, but the normalization uncertainty is quite large. Similar conclusions can be drawn for ϒ production, as can be seen from Fig. 2 where we compare our results with ALICE [22] and LHCb [39] data, although we observe that the uncertainty band of our calculation is narrower than for J=ψ production.
In Fig. 3 we show the comparison of our results for the J=ψ P ⊥ spectrum with ALICE and LHCb data. Again the agreement with the data is quite good except at very large transverse momenta where the CGC description is not expected to be very accurate. For ϒ production, as shown in Fig. 4, the agreement is also quite good at higher P ⊥ , but this leading order framework does not describe the increase in the mean P ⊥ of the ϒ compared to the J=ψ.

B. Proton-lead collisions
In Fig. 5 we show the differential J=ψ cross section as a function of the rapidity in pA collisions at a center of mass energy of 5 TeV and compare with ALICE [23] and LHCb FIG. 1 (color online). Differential J=ψ cross section as a function of Y in pp collisions. Data from Refs. [19,22]. [24] data. As in the pp case, the slope of the data is compatible with the one predicted by our calculation but the normalization uncertainty is large. In Fig. 6 we show the same quantity for ϒ meson production and compare our results with data from the ALICE Collaboration [40]. Here we can see that the normalization uncertainty of our calculation is again smaller than for J=ψ production, but there are only two data points to compare with, which limits the scope of this comparison.
In Fig. 7 we show the differential cross section as a function of the transverse momentum of the J=ψ in pA collisions at a center of mass energy of 5 TeV and compare it with LHCb [24] data. The agreement with the data is good but it is not possible to see if our calculation overestimates the cross section at large P ⊥ like it does for pp collisions since there is only one experimental bin between 8.5 and 14 GeV.
The normalization uncertainty from the charm quark mass and the factorization scale is quite large in the cross section itself. The modification of the production FIG. 3 (color online). Differential J=ψ cross section as a function of P ⊥ in pp collisions (2.5 < Y < 4). Data from Refs. [19,22].
FIG. 4 (color online). Differential ϒ cross section as a function of P ⊥ in pp collisions (2.5 < Y < 4). Data from Refs. [22,39].  2 (color online). Differential ϒ cross section as a function of Y in pp collisions. Data from Refs. [22,39]. mechanism due to the stronger color fields in nuclei is, however, very little affected by these factors. Indeed, if the nuclear geometry is properly treated within the Glauber picture, the CGC framework provides a very robust prediction for the nuclear modification factor, defined as This ratio has been extracted from the proton-proton and proton-nucleus cross sections by both the ALICE and the LHCb collaborations. In Fig. 8 we show it in the case of J=ψ production as a function of the rapidity and compare with experimental data. We observe that our results are close to the data when taking into account all sources of uncertainty. In particular, as already discussed in the Introduction, our calculation is much closer to the data than early CGC calculations in Ref. [25]. The nuclear modification factor for ϒ production was also measured by the ALICE [40] and LHCb [41] collaborations, in two and one rapidity bins respectively. We show the comparison of these data with our results on Fig. 9. Here no firm conclusions can be drawn at the moment since experimental uncertainties are quite large.
Finally we look at the behavior of R pA as a function of P ⊥ . Figure 10 shows the behavior of R pA as a function of P ⊥ for J=ψ production, compared to the measurement by the ALICE Collaboration [42]. Except at very low P ⊥ , the agreement of our calculation with this measurement is quite good, and the normalization uncertainty has almost completely canceled out. In Fig. 11 we also compare our calculation to PHENIX results [43] at lower energies. While our calculation displays a quite strong "Cronin" peak that is not seen in the data, also here our results are closer to data than the previous CGC calculations [25]. For the lower energy the calculation is more sensitive to the precise functional form of the initial condition. A pure MV model initial condition, which does not describe HERA data equally well, has a smaller Cronin enhancement. This feature was also noticed in Ref. [25]. The nuclear modification factor results at LHC energies, Figs. 8 and 10, are the main result of this paper.

V. DISCUSSION AND OUTLOOK
At high enough transverse momentum, particle production in QCD processes should be well described by collinear factorization. Indeed many comparisons of NLO calculations using nuclear parton distributions have been performed in the literature, see e.g. Ref. [44]. In general it seems that the shadowing present in the standard EPS09 [45] nuclear gluon distribution tends to underestimate the nuclear suppression observed experimentally in proton-nucleus collisions at lower transverse momenta. The collinear NLO calculation is in fact not that far from the physical picture of our work here. A leading order collinear calculation of heavy quark pair production results in a delta-function P ⊥ -distribution for the pair. To get a physically meaningful result for the differential cross section, one must go to NLO, where an additional radiated gluon can carry the recoil transverse momentum. In a standard NLO calculation, one integrates over the phase space of this gluon. At high energy this additional phase space integration can contribute a large logarithm ∼ ln s. The BK equation is, in this context, a resummation of these logarithmically enhanced terms to all orders in α s ln s. The large nuclear suppression in the data could be taken as an indication that the nonlinear effects included in BK indeed are significant for the J=ψ.
In "cold nuclear matter" (CNM) energy loss calculations such as Ref. [46], the physical argument behind the calculation can also be described in a language that is very similar to the one used here. The picture in the target rest frame is that of an incoming gluon splitting into a quark-antiquark pair. The pair can then lose energy as it propagates through the target nucleus, leading to a nuclear suppression of the cross section. This picture does not rely on the eikonal approximation, making it possible to extend the formalism to higher x A and transverse momenta. The practical implementation of the CNM model relies on a factorization of the quark pair production cross section and subsequent propagation through the nucleus. This factorization becomes questionable at very high energies in the target rest frame (very high s and forward rapidities), when the typical target transverse momentum scale Q s is of the order of the charm quark mass. In this limit the time scales of the quark pair formation and of the gluon radiation leading to the energy loss are not well separated. Indeed, in the CGC calculation the two are not factorized, but the pair production and propagation through the target nucleus are described coherently in the same formalism.
There is a normalization uncertainty in our calculation from the fact that we are using the dipole cross section [8] that was obtained for light flavors only. Requiring separately a good description of the charm DIS cross section F c 2 would require additional parameters in the fit, e.g. as a separate σ 0 for heavy quarks as in the AAMQS [2] parametrization. This would improve our confidence in the normalization of the charm production cross section, but would be a significant work outside the scope of this paper.
We also used the very crude color evaporation hadronization model to obtain the J=ψ cross section from the quark pairs. Our estimate for the normalization uncertainty was rather conservative, as we varied the quark mass without adjusting the value of the CEM constant F J=ψ together with it. The description of the transition from a quark pair to a J=ψ bound state has seen many recent theoretical advances in terms of nonrelativistic QCD (NRQCD) [47,48]. This has recently been successfully combined [49][50][51] with the small-x formulation for quark pair production. It would be interesting to combine the newer proton and nucleus dipole cross sections used here with the NRQCD hadronization picture. This is, however, left for future work.
Both of the uncertainties discussed above affect more seriously the separate proton-proton and proton-nucleus cross sections. They should mostly cancel in the nuclear modification factor R pA . We leave a more detailed investigation of the separate cross sections to future work. The main purpose of this paper has been to focus on the importance of a proper treatment of the nuclear geometry, including the use of a proton transverse area consistent with the DIS fits used to obtain the dipole cross section. We have shown that when this is done using the dipole cross section parametrization obtained in Ref. [8], the resulting nuclear modification ratio is much closer to the experimental data than in previous CGC estimates.  11 (color online). Nuclear modification factor for J=ψ production at RHIC energies as a function of P ⊥ at Y ¼ 1.8. PHENIX data (1.2 < Y < 2.2) from Ref. [43].