Heavy-quarkonium suppression in p-A collisions from parton energy loss in cold QCD matter

The effects of parton energy loss in cold nuclear matter on heavy-quarkonium suppression in p-A collisions are studied. It is shown from first principles that at large quarkonium energy E and small production angle in the nucleus rest frame, the medium-induced energy loss scales as E. Using this result, a phenomenological model depending on a single free parameter is able to reproduce J/psi and Upsilon suppression data in a broad xF-range and at various center-of-mass energies. These results strongly support energy loss as the dominant effect in heavy-quarkonium suppression in p-A collisions. Predictions for J/psi and Upsilon suppression in p-Pb collisions at the LHC are made. It is argued that parton energy loss scaling as E should generally apply to hadron production in p-A collisions, such as light hadron or open charm production.


Introduction
The spectacular quenching of hadrons produced at large p ⊥ in Pb-Pb collisions at the LHC [1,2], as well as the jet imbalance reported in those collisions [3,4], find a natural explanation in terms of parton energy loss in a quark-gluon plasma (QGP). For light hadron production at mid-rapidity and sufficiently large p ⊥ , the parton energy loss is dominantly radiative, in average of the form ∆E ∼ α sqhot L 2 [5,6], with L the distance travelled by the parton through the hot medium andq hot the rate per unit length of transverse momentum broadening in the medium. The strength of jet-quenching can be explained if the transport coefficientq in a QGP is larger, by one to two orders of magnitude, than its estimate in cold nuclear matter,q cold ∼ 0.045 GeV 2 /fm [7]. This is why jet-quenching is considered as a prominent QGP signal. However, despite the wealth of data accumulated so far at RHIC and LHC, the in-depth understanding of energy loss processes in a QGP remains far from complete (see [8] for a discussion).
Drastic nuclear suppression effects are not only seen in A-A but also in p-A collisions, at least for some processes and in some kinematical conditions. For instance, quarkonium [9] but also light hadron [10,11] production at large longitudinal momentum fraction x F (or large rapidity) is strongly suppressed in p-A as compared to p-p collisions. Understanding nuclear suppression in cold nuclear matter, a well-controlled medium as opposed to an expanding QGP, should be a prerequisite in order to interpret quantitatively nuclear suppression in heavy-ion collisions. However, it is striking that there is no consensus yet on the origin of J/ψ suppression at large rapidity/x F in p-A collisions, from SPS to RHIC [9,12,13], despite many theoretical attempts (see [14] for a review).
Recently, new scaling properties have been identified for the induced gluon radiation spectrum dI/dω, and associated energy loss ∆E, of hard processes where a color charge undergoes small angle scattering through a static medium (cold matter or QGP) [15]. In the present work we address the phenomenological consequences of these results on J/ψ and Υ nuclear suppression in p-A and π-A collisions, parametrized by the attenuation factor (in the following we use the generic notations "ψ" and "p-A") We will show that the available large-x F quarkonium suppression data in p-A collisions can be explained by parton energy loss in cold matter. Althoughq in cold matter is small, a strong nuclear attenuation arises due to the specific parametric behaviour ∆E ∝ E at sufficiently large E, where E is the quarkonium energy in the target nucleus rest frame. As discussed in Ref. [15] and reviewed in the present paper (Section 2), this behaviour holds when the hard partonic subprocess can be viewed, in the nucleus rest frame, as the small angle scattering of a color charge. In the following we focus on quarkonium hadroproduction (see Fig. 1a), where the heavy quark mass provides the hard scale allowing for a perturbative QCD description, and for which p-A suppression data are quite abundant. Our discussion should however apply more generally to hadron hadroproduction, for instance to light hadron production in p-A collisions (provided the light hadron p ⊥ plays the role of the hard scale, i.e., p ⊥ ∼ > 1 GeV), see Fig. 1b. Light hadron nuclear suppression due to parton energy loss will be addressed in a future work.
As is well-known, the quarkonium hadroproduction mechanism in elementary p-p collisions is still under debate 1 (see for instance Ref. [18] for a review). In order to study (c) (d) quarkonium nuclear suppression in the most model-independent way, in the present study we will only assume that the heavy-quark QQ pair of mass M is produced, within the perturbative proper time scale τ QQ ∼ 1/M , in a compact color octet state, and remains color octet for a time τ octet ≫ τ QQ . In quarkonium production models where color neutralization is a soft, non-perturbative process, τ octet coincides with the quarkonium hadronization time τ ψ , and this assumption holds at any x F . In the Color Singlet Model (CSM), the gluon emission required for color neutralization of the QQ pair is constrained at large enough x F to become softish and thus to occur late, leading to τ QQ ≪ τ octet τ ψ . Thus, the assumption of a color octet QQ pair living longer than the perturbative time scale ∼ 1/M holds quite independently of the quarkonium production model. 2 Working in the nucleus rest frame and considering the limit E ≫ M ∼ > p ⊥ , quarkonium hadroproduction looks like small angle scattering of an "asymptotic" color charge, i.e., prepared in the "far past" and propagating in the "far future" as compared to the perturbative time scale t hard = τ QQ · (E/M ) ∼ E/M 2 . This is illustrated in Fig. 1a for the generic gg → QQ partonic subprocess, viewed in the nucleus rest frame as the splitting g → QQ of the incoming gluon, followed by a rescattering in the target. 3 In the present study we will assume for simplicity that the octet QQ pair arises dominantly from the splitting of an incoming gluon. This should be a valid assumption for all p-A data considered in this paper, except at very large values of x F (x F ∼ > 0.8), where quark-induced processes come into play. (We will further comment on this point in Section 4.1.) Within those assumptions, the associated gluon radiation with large formation time, t f ≫ t hard , is similar to (non-abelian) Bethe-Heitler radiation off a fast color octet undergoing an effective transverse momentum kick q ⊥ . The typical q ⊥ is expected to be larger in p-A than in p-p collisions due to transverse momentum broadening being proportional to the target size L, ∆q 2 ⊥ ≡ ℓ 2 ⊥ ≃q L. 4 As a result, the medium-induced radiation spectrum is similar to the Bethe-Heitler spectrum (2.8), up to the replacement of the total momentum transfer by the broadening ℓ 2 ⊥A through the nucleus A (see (2.17) and Section 2), with M ⊥ = (M 2 + p 2 ⊥ ) 1 2 the transverse mass of the QQ pair and Λ 2 p = max(Λ 2 QCD , ℓ 2 ⊥p ). This leads to an average loss ∆E ∝ E. When Λ 2 p < ℓ 2 The spectrum (1.2) is at the basis of the phenomenological study presented here, whose main results can already be found in Ref. [19]. The scaling ∆E ∝ E in quarkonium hadroproduction was first postulated in [20] (also revisited in [21]) yet this assumption was not motivated and the parametric dependence on L and M arbitrary (and different from (1.3)). In Ref. [22], an energy-independent bound on ∆E was derived, but in a specific setup where the nuclear broadening of the final tagged particle was neglected (we will briefly comment on this point in Section 2).
We stress that the spectrum (1.2) is coherent. This can easily be seen in a calculation using physical polarizations for the radiated gluon (see Section 2). With this choice the medium-induced radiation spectrum indeed arises from the interference between the initial and final state emission amplitudes. The energy loss (1.3) is thus neither a purely initial nor final state effect, and is distinct from gluon radiation resummed in leading-twist parton distribution and fragmentation functions. Being process-dependent (e.g., it is expected in J/ψ production in p-A collisions but not in Drell-Yan production, see below) and suppressed by a power of the hard scale M ⊥ , it is naturally interpreted as a higher-twist effect. Nevertheless it plays a crucial role in a broad x F or rapidity interval, as we shall see.
To make the physics under consideration clear, let us mention that the spectrum (1.2) is not expected in quarkonium (real or virtual) photoproduction, nor in inclusive deep inelastic scattering (DIS) off nuclei, where the incoming energetic particle participating to the hard subprocess is colorless (see Fig. 1c for the DIS case). 5 In those cases radiation with t f ≫ t hard can only arise from final state radiation. The latter (DGLAP-like) radiation is independent of the medium properties and cancels in the medium-induced 4 In our study the broadening ∆q 2 ⊥ equals the amount of soft rescattering ℓ 2 ⊥ . In other words it is defined with respect to an "ideal" target where soft rescatterings are absent, see Section 2.2.
5 This is to be distinguished from resolved photoproduction, which at the partonic level is similar to hadroproduction (Fig. 1a), and where we thus expect ∆E ∝ E. spectrum. Quite remarkably, at large z the quarkonium production data in deep inelastic muon scattering [23] exhibit no nuclear suppression (but instead a slight enhancement), in sharp contrast to hadroproduction.
Drell-Yan (DY) production off nuclei is similar to quarkonium photoproduction, since by definition the energetic particle produced perturbatively (the Drell-Yan photon of mass Q) is colorless, see Fig. 1d. Radiation with t f ≫ t hard ∼ E/Q 2 must arise from initial state and does not contribute to the medium-induced spectrum. Thus, neither the spectrum (1.2), nor the energy loss ∆E ∝ E of the type (1.3), is expected in DY production. This is qualitatively consistent with the much milder nuclear suppression of DY production [24] when compared to J/ψ hadroproduction in the same kinematical range.
In Section 2 we discuss the medium-induced gluon radiation spectrum, derived in [15], which is at the basis of the model detailed in Section 3. Phenomenological applications are presented in Section 4 and we conclude the paper by a discussion (Section 5).

Revisiting energy loss scaling properties
In this section we justify the expression (1.2) for the gluon radiation spectrum associated to quarkonium hadroproduction, in a less heuristic way than in Ref. [15].

Asymptotic charge
We first consider the case of an on-shell ("asymptotic") parton of energy E undergoing an elastic scattering and exchanging a gluon with transverse momentum ℓ ⊥ with a nuclear target, see Fig. 2a. As is well-known, scattering can induce radiation, provided the quantum state of the charge is perturbed. In QED, this happens when the scattering angle θ s ≃ ℓ ⊥ /E is non-zero. In QCD, radiation can occur even in the limit θ s → 0, due to the incoming parton color rotation in the elastic scattering. These features are illustrated by the expression of the gluon radiation amplitude induced by the elastic scattering, given by the diagrams of Fig. 2b. We denote by ω and k ⊥ the radiated gluon energy and transverse momentum. We focus on soft (ω ≪ E) and small angle (k ⊥ ≪ ω) gluon radiation. For an on-shell light quark the radiation amplitude where θ M ≡ M/E and the lower cutoff arises from the constraint k ⊥ > Λ (with Λ ≡ Λ QCD ), put by hand for the consistency of our perturbative QCD treatment (we also assumed ℓ ⊥ > Λ). We obtain the soft radiation spectrum off an on-shell quark of mass M [15,26], The average energy loss of a heavy quark,

Color charge resolved in a hard process
The case of a color charge resolved in a hard process can be simply illustrated by inserting a hard exchange q ⊥ ≫ ℓ ⊥ in the light quark scattering process, see Fig. 3a. The transfer ℓ ⊥ now plays the role of nuclear momentum broadening, ℓ 2 ⊥ ≡ ∆q 2 ⊥ . Since some radiation is released even when ℓ ⊥ = 0 due to the presence of the hard exchange, the relevant quantity is the medium-induced radiation spectrum, We concentrate on soft radiation as compared to the hard process, ω ≪ E, k ⊥ ≪ q ⊥ , and on the limit of large gluon formation time. It can indeed be checked that the domain t f ≫ max(L, t hard ) contributes most to the medium-induced loss (see Section 2.3). The radiation amplitude is thus dominated by the diagrams of Fig. 3b. Diagrams of the type shown in Fig. 4 can be neglected. The diagram of Fig. 4a, relevant when t f L, contributes to the energy loss as ∆E ∝ L 2 , which is negligible compared to the contribution ∆E ∝ E we focus on throughout our study. The diagram of Fig. 4b is suppressed due to ℓ ⊥ ≪ q ⊥ .
Note that in the limit t f ≫ L, the induced radiation cannot probe the nuclear size L, but is however sensitive to the total amount of soft rescattering, ℓ 2 ⊥ =qL, transferred to the fast color charge, independently of the actual number of soft rescatterings. Modeling soft rescatterings by a single effective (semi-hard) scattering with ℓ 2 ⊥ =qL should thus allow addressing the main features of medium-induced radiation. The azimuthal integral yields [15] where x ≡ ω/E and we added the constraint k ⊥ > Λ in the k ⊥ -integral. Approximating It is easy to check that for the scattering of a fast octet charge (the case of interest), the color factor F c in (2.15) reads F c = N c . Note that the virtual correction shown in Fig. 5b contributes a factor 2 in (2.15). 9 We can also verify that with a parton mass M = 0 we get the same expression as (2.15), up to the change in the hard scale q ⊥ → M ⊥ = (M 2 + q 2 ⊥ ) 1 2 [15]. The spectrum off an energetic color octet charge thus reads where the Θ-function reminds us that only ℓ ⊥ > Λ can induce the emission of perturbative gluons. Quite remarkably, the medium-induced spectrum (2.16) is parametrically similar to the radiation spectrum (2.8) of an asymptotic parton of "mass" M ⊥ . Finally, we stress that in the above calculation, the medium-induced spectrum has been defined with respect to an "ideal" target for which ℓ ⊥ = 0, see (2.10). In practice, nuclear suppression is measured in a nucleus A with respect to that in a nucleus B < A, where soft rescatterings can also occur (even in the proton case B = 1). The mediuminduced spectrum relevant to this situation should be defined with respect to the nucleus B and reads (as already announced in (1.2) for B = 1) where Λ 2 B ≡ max(Λ 2 QCD , ℓ 2 ⊥B ). When the broadening in the nucleus B is too small, ℓ 2 ⊥B < Λ 2 QCD , the spectrum (2.17) coincides with the spectrum (2.16) defined with respect to the "ideal" target with ℓ ⊥ = 0. When ℓ 2 ⊥B > Λ 2 QCD , the induced radiation in A with respect to B becomes independent of Λ QCD .

Application to quarkonium hadroproduction
As illustrated by the previous section, the behaviour ∆E ∝ E for medium-induced parton energy loss is not forbidden by first principles. We expect such a behaviour in all p-A processes where some color charge is scattered at small angle in the nucleus rest fame, in particular in quarkonium hadroproduction at large E.
However, in order to apply the spectrum (1.2) to quarkonium hadroproduction, the underlying partonic subprocess should effectively look like the scattering of an energetic, pointlike color charge, at least within the formation time t f ≫ L of the radiated gluon. This is the case if where t hard is the hard process time scale, t octet the lifetime of the color octet QQ pair, and t ψ the quarkonium hadronization time (see the Introduction). The second condition states that the QQ pair is effectively pointlike when the transverse wavelength 1/k ⊥ of the radiated gluon is larger than the transverse size r ⊥ of the quark pair at a time ∼ t f . In principle, the two conditions (2.18) can be checked a posteriori, using the typical formation time contributing to the observable of interest (in our case nuclear attenuation).
As an illustration, we show that the typical t f contributing to the average loss (1.3) formally satisfies (2.18). For simplicity we assume here L < t hard ∼ E/M 2 ⊥ and M ⊥ ≃ M , and denote ℓ 2 ⊥ ∼ ∆q 2 ⊥ . It is trivial to check that ∆E in (1.3) arises from radiated energies ω ∼ E (ℓ ⊥ /M ) ≪ E and transverse momenta k 2 ⊥ ∼ ℓ 2 ⊥ (see (2.14)). The typical t f thus satisfies the inequalities arising from the nuclear broadening ℓ ⊥ being soft compared to M , and hard (for large nuclei) compared to the non-perturbative scale τ −1 ψ ≃ 0.6 GeV (see Section 3.5.1). In this respect, let us recall that momentum broadening is related to the saturation scale Q s in the nucleus, ℓ 2 ⊥ = Q 2 s [28], indeed considered as a semi-hard scale. The second condition in (2.18) reads where the relative transverse velocity of the heavy quarks is estimated by v ⊥ ∼ p B⊥ /E, with p B ∼ α s M the Bohr momentum of the quarkonium state. 10 Thus, the conditions (2.18) are fulfilled in the perturbative domain α s ≪ 1 and provided nuclear broadening is a semi-hard scale. This defines the theoretical limit where the energy loss (1.3) can be applied to quarkonium hadroproduction. In practice, ℓ ⊥ ∼ √q L ≫ τ −1 ψ might not be satisfied. Indeed, usingq = 0.08 GeV 2 /fm (see Section 4) and L = 7 fm we have √q L ≃ 0.7 GeV ∼ τ −1 ψ , somewhat questioning the validity of (1.3) for quarkonium hadroproduction. However, it should be stressed that the observable of interest in this paper -nuclear attenuation -does not directly depend on the average energy loss, but rather on the energy loss probability distribution P(ω, E), see Section 3. As is well-known and generic to jet-quenching phenomenology [7], nuclear attenuation is dominated by the low energy tail of P(ω, E), i.e., by ω much smaller than the typical ω contributing to the average loss. This leads to smaller values of t f to be used in (2.19), leading to the required condition t f ≪ t ψ .
Under the conditions (2.18) we thus expect the induced radiation spectrum to be given by Eq. (1.2), as derived in Section 2.2 for a fast color octet charge. It should be clear from Section 2.2 that the parametric dependence of the spectrum is uniquely determined, and should thus apply to other processes than quarkonium production (like open charm and light hadron hadroproduction), as discussed in the Introduction.
Finally, we note that within our approximation (2.18) the quarkonium bound state is formed far beyond the nucleus, t ψ ≫ L. This approximation may break down at low proton beam energy or at small values of x F . We will comment on this when discussing the limits of applicability of the model in Section 3.5.1.

Shift in energy or "medium-induced splitting"
The starting point of the model consists in expressing the J/ψ differential production cross section dσ/dE in p-A collisions simply as that in p-p collisions, with a shift in the quarkonium energy E accounting for the energy loss ε incurred by the octet cc pair propagating through the nucleus, The energy loss ε is more conveniently defined in the nucleus rest frame, and we thus denote E and E p ≃ s/(2m p ) the J/ψ and projectile proton energies in this frame (with m p the proton mass). We have ε ≤ E p − E from energy conservation, and we impose ε ≤ E for consistency with the soft radiation approximation. Hence ε max = min(E p − E, E) in (3.1). The quantity P(ε, E) is the energy loss probability distribution or quenching weight associated to the radiation spectrum (1.2), to be discussed in Section 3.4. The measured differential cross sections are usually expressed as a function of x F (or of the rapidity) rather than E. The variable x F is defined as the longitudinal momentum fraction between the J/ψ and projectile proton in the c.m. frame of an elementary p-N collision (of energy √ s). It can be related to the J/ψ transverse mass M ⊥ and rapidity y ′ in this frame. In the limit √ s ≫ m p , where cosh ∆y = √ s/(2m p ) with ∆y the projectile proton rapidity in the c.m. frame, we obtain from (3.2) The relation (3.3) can be inverted to give Introducing the variable x ′ , in between brackets in (3.6) approaches unity, and the p-p cross section is evaluated at a shifted value of is thus equivalent to a simple translation in x F . This is not true at all values of x F , due to the presence of the Jacobian. In the following we will use the expression (3.6), where E = E(x F ), and the relations (3.3) and (3.4), valid at all x F . In a model where the J/ψ is produced through a 2 → 1 partonic subprocess, the expression (3.3), denoted as E ≡ x 1 E p , simply arises from the standard relations between parton momentum fractions, It is interesting to mention that the main equation of our model (3.1) is equivalent to where z ≡ E/(E + ε) is interpreted as a (medium-induced) splitting variable describing the energy loss process (E is the energy of the charge after radiating the energy ε), and F loss (z) as a "medium-induced splitting function". The expression (3.7) follows from (3.1) by changing variable from ε to z (giving z min = max(E/E p , 1/2)), and using the fact that the quenching weight P(ε, E) is a scaling function of the ratio ε/E, 11 In writing (3.1), we implicitly assumed that the energy of the radiating octet cc pair is the same as the final J/ψ energy E. However, the equivalent expression (3.7) suggests that our nuclear suppression model based on a simple shift in E might apply to more general situations where the final detected particle's energy arises from the fragmentation of some parent parton's energy, with a fragmentation variable z < 1. Indeed, suppose that the observable p-p cross section is of the form 12 and that medium-induced radiation and hadronization factorize, 13 where z ′ min = max(E/(zE p ), 1/2). Changing the order of the z and z ′ integrals in (3.10) and using (3.9), we recover (3.7).
The quarkonium p-A cross section is thus related to the observable p-p cross section according to (3.7), or equivalently (3.1), independently of the form of the fragmentation function D(z). This result mostly follows from the scaling property (3.8) of the quenching weight, and is thus expected to hold in all processes involving a fractional medium-induced energy loss (∆E ∝ E), in particular in open charm and light hadron hadroproduction. Instead of describing the energy loss process as a shift in energy (see (3.1) or (3.6)), one could alternatively describe it as in (3.7) by a medium-induced splitting function 14 F loss (z).
Finally, although in the above discussion we assumed the p-p cross section to obey factorization (see (3.9)), we believe that (3.1) might hold independently of this assumption. For instance, as long as the underlying partonic process is similar to the scattering of a color charge (as in Figs. 1a and 1b), we may imagine the p-p cross section at large E to be affected by late comover rescattering [30] and thus to violate factorization, and nevertheless the p-A cross section to be given by (3.1). The only crucial assumption is that the partonic subprocess induces radiation, as dictated by perturbative QCD, independently of the precise mechanism fixing the quantum numbers of the final detected particle.

Absolute production cross section
The dynamics of heavy-quarkonium production in hadronic collisions is still uncertain. In particular, none of the existing models proposed to describe heavy-quarkonium production is able to reproduce simultaneously all the features reported experimentally, at both p ⊥ M and p ⊥ ≫ M .
In the present approach, a crucial ingredient entering (3.6) is the x F single differential absolute cross section, dσ ψ pp /dx F , of J/ψ and Υ production in p-p collisions at a given center-of-mass energy. In order to be as model-independent as possible, dσ ψ pp /dx F used in (3.6) is not taken from theory but determined from a fit to the data. We found that it can be conveniently parametrised as where x ′ is defined in (3.5) and the exponent n is extracted from p-p 15 and π − -p data taken at the same center-of-mass energy -whenever possible -as the p-A and π − -A 13 This should be guaranteed by the separation of time scales, t f ≪ t ψ , see Section 2.3. 14 This designation is motivated by the fact that F loss is perturbatively calculable and can be Mellin convoluted with the (vacuum) fragmentation function D(z) to give the "medium modified fragmentation function" D med (z) = dz ′ D(z ′ ) F loss (z/z ′ ). We however stress that F loss (z) is process-dependent (e.g., it is present in quarkonium but absent in DY production, see the Introduction). Thus, D med (z) in the latter equation differs from the medium modified fragmentation functions assumed to be universal and discussed elsewhere, see for instance Ref. [29]. 15 Table 1. Values of n extracted from J/ψ production in p-p (or p-Be, p-C) and π − -p collisions; see text for details. Table 2. Values of n extracted from Υ production in p-p (or p-d) collisions; see text for details. measurements discussed in this paper. 16 Note that the normalization parameter in (3.11), or equivalently the total production cross section, is irrelevant for our purpose since only nuclear production ratios are considered in this paper, see Eq. (1.1). The values of n for J/ψ production from SPS to LHC energy are summarized in Table 1. 17 The index for p-p production grows smoothly from n ≃ 4-5 at low energies ( √ s 40 GeV) up to n ≃ 8 at RHIC and n ≃ 30 at the LHC. At LHC, the relative uncertainty on n is as large as δn/n ∼ 20% (n = 32.3 ± 7.5) because of the too small x F domain, |x F | 0.02, covered by the data and the fact that around mid-rapidity, x ′ ≪ 1, the parametrization (3.11) becomes independent of n. However we checked that the resulting uncertainty on R J/ψ pA at LHC (and similarly at RHIC) is marginal because of these very reasons. In the Υ channel, the exponents are slightly smaller than for J/ψ production; see Table 2.
The comparison between the fits and J/ψ (respectively, Υ) production data is shown in Fig. 6 (respectively, Fig. 7). An excellent agreement is observed on a very wide range of center-of-mass energies, spanning from √ s = 16.8 GeV to √ s = 7 TeV, on the full x F domain covered experimentally, and for both J/ψ and Υ production. These results indicate that the parametrization (3.11) can be safely used in (3.6) in order to compute heavy-quarkonium nuclear suppression. Independently of the present work, this parametrization can prove useful in future phenomenological studies for which a data-driven knowledge of the x F single differential cross sections is necessary.  Table 1. Data are taken from [12,[31][32][33][34].  Table 2.

Transport coefficient and nuclear broadening
The amount of medium-induced gluon radiation, and hence the strength of ψ suppression in p-A collisions, is controlled by the nuclear broadening ∆q 2 ⊥ ≡ ℓ 2 ⊥ in Eq. (1.2). For a path length L crossed in the target, the broadening reads where the transport coefficientq in the target nucleus is related to the gluon distribution G(x) in a target nucleon as [38] 18 In the latter equation ρ is the target nucleon number density, and the scaling violations in the running of α s and in the evolution of the gluon density are neglected sinceqL 1 GeV 2 .
The typical value of x at which xG(x) should be evaluated in the r.h.s. of (3.13) is discussed in Appendix A. The result depends on whether the hard subprocess is incoherent, t hard ≪ L, or coherent, t hard ≫ L, in the nucleus. Assuming 2 → 1 subprocess kinematics , so that the incoherent and coherent limits correspond respectively to Using the power-law behavior xG(x) ∼ x −0.3 suggested by small-x (x < 10 −2 ) fits to HERA data [39], we can thus extract the x and ρ dependence ofq by writinĝ where ρ 0 is in principle an arbitrary constant density, andq 0 ≡q(x = 10 −2 , ρ = ρ 0 ). In Ref. [19] we used the hard sphere (HS) approximation 19 and thus the same uniform density ρ = ρ HS = [(4/3)πr 3 0 ] −1 ≃ 0.17 fm −3 for all nuclei. Within this approximation, q =q(x) withq 0 ≡q(x = 10 −2 , ρ = ρ HS ), the average of L = dz is found to be L HS = 3R A /2, and the broadening is directly obtained from (3.12). In the present study we will use more realistic (non-uniform) nuclear density profiles ρ(r) extracted from electronproton scattering experiments [40]. In order to formally recover the situation considered in Ref. [19] when ρ(r) is constant, we thus choose ρ 0 = ρ HS in (3.14).
When ρ(r) is not constant, the parameter L = dz entering (3.12) is badly defined, but the broadening is still well-defined, since it is proportional to ρ dz rather than to dz. Using (3.14) we can write The effective path length L eff is mathematically well-defined and can be related to the number N part of nucleons participating to the broadening of the fast color charge. Using dN part = ρ σ dz, where σ is interpreted as the cross section for having non-zero broadening in parton-nucleon scattering, we obtain 20 For minimum bias p-A collisions, the average of N part in the events with J/ψ production can be calculated within Glauber theory and reads where we used the normalization d 3 r ρ(r) = d 2 b T A (b) = A. The effective path length becomes  Table 3. Values of L eff for various nuclei, using (3.19) and realistic nuclear densities [40], or within the hard sphere approximation.
As can be seen from Eq. (3.18), the additional (effective) path length in a nucleus with respect to a proton is independent of σ, and can thus be uniquely determined knowing the nuclear density profile. For the effective length in the proton, we take L p = 3R p /2 = 1.5 fm, using R p = 1 fm for a generic proton length scale. 21 In summary, L eff is given by The values of L eff obtained from (3.19) using realistic nuclear density profiles [40] are listed in Table 3 for various nuclei, and compared to the values L HS = 3R A /2 corresponding to the hard sphere approximation previously used in Ref. [19].

Energy loss probability distribution
In this section we discuss the quenching weight P(ε, E) entering (3.6). A well-known procedure (used for instance in the case of large p ⊥ jet-quenching in A-A collisions [7]) to construct a normalized P(ε, E) from the single gluon emission spectrum dI/dω, consists in assuming independent emissions of soft gluons. In this so-called Poisson approximation, Let us mention that P(ε, E) is a solution of the equation where by convention P(ε < 0, E) = 0, and L stands for any parameter entering the expression of dI/dω. If L is the medium length crossed by the fast charge, (3.21) is formally identical to the kinetic equation used by Landau to study ionization losses in normal matter [41]. An important feature of the Poisson approximation (3.20) is that not only each ω i , but also the accumulated loss ω i , is supposed soft as compared to E. In other words, the "energy degradation" of the fast particle during the multiple emission process is neglected. Whether the Poisson approximation is appropriate or not obviously depends on each energy loss process and on the specific properties of dI/dω, in particular on the multiplicity N ω of radiated gluons with energy ∼ O (ω). In the present case, dI/dω given in (2.17) is a scaling function of ω/ω A , whereω A ≡ (ℓ ⊥A /M ⊥ ) E, 22 and N ω is estimated as When ω ≪ω A , N ω becomes potentially large, α s ln (ω A /ω) ∼ O (1), questioning the assumption of independent multiple emissions. 23 When ω ∼ >ω A , N ω is small, N ω O (α s ) ≪ 1. However in this region each emitted gluon carries away a fixed fraction ∼ ℓ ⊥A /M ⊥ of the energy E, and the fast particle energy degradation cannot be neglected. 24 Thus, in our context the Poisson approximation proves fishy.
A simple way to deal with this problem is to supplement (3.20) with the condition that the energy ε is carried away by a single gluon, i.e., δ (ε − ω i ) → n δ (ε − ω j ) in (3.20). This yields the (normalized) quenching weight The latter is simply interpreted as the product between the "probability" dI/dε to radiate a gluon with ω j = ε and the probability (given by the exponential Sudakov factor) to have no extra radiation carrying ω k ∼ > ε. In our context, the expression (3.23) of the quenching weight is better founded than the Poisson expression (3.20), and will be used in (3.6).

Other nuclear effects
Besides energy loss effects, other mechanisms might affect ψ suppression in nuclei. In this section, the role of nuclear absorption and saturation on ψ suppression in p-A collisions is discussed.

Nuclear absorption
At small ψ energy in the nucleus rest frame, the hadronization time t ψ = τ ψ · (E/M ) (where τ ψ is the proper hadronization time) becomes comparable to the typical nuclear size, t ψ L. Consequently, ψ states are produced on average within the target nucleus and might suffer inelastic interaction with nuclear matter, the so-called nuclear absorption process. From (3.3), this should be the case at low proton beam energy E p (i.e. at low √ s ≃ 2m p E p ) or at small values of x F . The J/ψ suppression in p-A collisions at the 22 For the present discussion we neglect the second term of the spectrum (2.17). 23 In the case of large p ⊥ jet-quenching [7], a large multiplicity at small ω is compensated by a factor t f /L ≪ 1, resulting in a small gluon occupation number ∼ (t f /L) Nω ≪ 1, thus supporting the assumption of independent emissions. In our case the spectrum (2.17) arises from large formation times t f ≫ L, and the estimates of gluon occupation number and gluon multiplicity coincide. 24 This problem was previously addressed in Ref. [42] in the context of large p ⊥ jet-quenching.
SPS (e.g. by the NA60 experiment at E p = 158 and 450 GeV [43]) has in particular often been attributed to nuclear absorption effects. Nuclear absorption effects are not included in this analysis for two reasons. First of all, the strength of nuclear absorption strongly depends on the (effective) absorption cross section σ ψ abs which is poorly constrained from data [44,45]. Moreover, when t ψ ∼ L, the hierarchy (2.18) upon which the medium-induced spectrum (2.17) relies is no longer valid and hence the use of the latter becomes dubious. In this paper we will therefore focus on the region t ψ ≫ L (the region of validity of (2.17)), where the effect of nuclear absorption is irrelevant.
In the calculations presented in Section 4, we shall indicate by an arrow the value of x crit F (or y crit ), defined as E(x crit F )/M × τ ψ = L, below which nuclear absorption might play a role. For the numerical values of x crit F , the J/ψ hadronization time is given by the mass splitting between 1S and 2S states, τ J/ψ = (M ψ ′ − M J/ψ ) −1 ≃ (0.6 GeV) −1 ≃ 0.3 fm (a similar estimate is obtained in the Υ channel). Note that x crit F becomes negative at large collision energy, in which case the model should not only apply at large positive x F but also down to x F < 0.

Saturation and nuclear PDF effects
At small values of x, partons inside the nucleus wavefunction start to overlap, leading to the phenomenon of saturation (see for instance [46] for a review). Although saturation effects should also occur in a proton, they are expected to scale roughly like the nucleus transverse density, V /S ∼ A 1/3 , therefore being stronger in large nuclei at a fixed value of x. As a consequence, the ψ normalized yield in p-A collisions is likely to be suppressed with respect to that in p-p collisions -independently of the energy loss effects discussed above -either at large x F and/or at high energies (RHIC, LHC) where small values of x are probed in the target nucleus.
The effects of (gluon) saturation on J/ψ suppression in p-A and A-A collisions have been addressed by many authors, see e.g. [47,48]. In the present paper, we shall implement the physics of saturation following the work of Fujii, Gelis and Venugopalan [48], where J/ψ suppression has been computed within the Color Glass Condensate assuming 2 → 1 kinematics for the production process. The nuclear suppression is a scaling function of the saturation scale Q s and can be simply parametrised as [48] S J/ψ with b = 2.65 GeV 2 and α = 0.417. Unfortunately, no equivalent parametrization has been given in the Υ channel. We will assume in the present approach that saturation effects on heavy-quarkonium production are a scaling function of Q s /M ⊥ . 25 Therefore, the Υ suppression due to saturation reads In order to make reliable predictions at RHIC and LHC, the J/ψ and Υ nuclear production ratio is determined assuming energy loss effects, R E.loss pA from Eq. (3.6), with and without saturation effects, The saturation scale appearing in (3.24) is closely related to the transport coefficient q given by (3.13), namely [28,49] In other words Q 2 s is nothing but the transverse momentum broadening discussed in Section 3.3, see Eqs. (3.12) and (3.15). The inclusion of saturation effects thus does not require any additional parameter once the parametrization (3.14) forq is employed andq 0 is determined.
Let us mention thatq 0 is related to the saturation momentum in a proton at x = 10 −2 , Q s0 . We have, from (3.14) and (3.26), withq 0 in GeV 2 /fm. Comparing the value ofq 0 obtained in our model from a fit to the E866 J/ψ nuclear suppression data and the current estimates of Q 2 s0 obtained from a fit to small-x 2 DIS data [50] should provide a non-trivial (though not conclusive) test of our model.
Another, earlier approach in order to model the modifications of parton densities in nuclei is the use of leading-twist nuclear PDF (nPDF) which have been determined from global fit analyses of e-A DIS or p-A Drell-Yan data for more than a decade (see e.g. [51]). In this framework, ψ production in p-A collisions is proportional to the gluon distribution in the nucleus G A (x 2 , M ⊥ ). Therefore, ψ suppression can be modelled as The predictions to be discussed in the next section will be performed assuming energy loss effects, supplemented with predictions including saturation effects at RHIC and LHC energies where these are expected to play a role. For completeness, we will also critically compare in Section 4.7 these results with those obtained using nPDF.

Phenomenology
After the description of the energy loss model, the phenomenology of ψ suppression in hadron-nucleus collisions is investigated in this section. In the practical applications, we take Λ QCD = 0.25 GeV, p ⊥ = 1 GeV in the transverse mass M ⊥ = M 2 + p 2 ⊥ , and M = 3 GeV (M = 9 GeV) for the mass of a compact cc (bb) pair. As we can easily verify a posteriori, the typical scale entering the running coupling constant is not too large, qL ∼ 1 GeV 2 , which justifies the assumption of a frozen coupling, α s = 1/2, at such semi-hard scales.

Fitting procedure
The only parameter of the model, the transport coefficientq 0 , is determined by fitting the J/ψ suppression measured by E866 [9] in p-W over p-Be collisions ( √ s = 38.7 GeV). This choice is motivated by the fact that the E866 measurements are the most precise performed so far and cover a wide range in x F . We choose to perform the fit in the [0.2-0.8] x Frange for the following reasons: at x F 0.2 J/ψ suppression might be affected by nuclear absorption (see Section 3.5.1) while at x F 0.8, we expect quark-induced subprocesses to come into play, possibly modifying the overall normalization of the medium-induced spectrum (1.2). Note also that this x F -range at E866 energy corresponds to values of x 2 10 −2 for which saturation effects are expected to be small, of the order of 5% at most on the W/Be ratio.
The fit givesq 0 = 0.075 ± 0.005 GeV 2 /fm, where the quoted uncertainty is determined from the χ 2 minimization procedure. A systematic uncertainty on the value ofq 0 can be roughly estimated by restricting the x F -range used for the fit to the interval [0.3-0.7]; we found that it would increase the value ofq 0 toq 0 ≃ 0.087 GeV 2 /fm.
The result of the fit is shown in Fig. 8 where excellent agreement is observed in the whole fit range. Note however that for x F 0.1, nuclear absorption is expected to play a role; see the vertical arrow at x crit F ≃ 0.07 below which the J/ψ formation time becomes smaller than the size of the target tungsten nucleus (see Section 3.5.1). It is worth mentioning that the fitted transport coefficient,q 0 = 0.07-0.09 GeV 2 /fm, would correspond to the saturation scale in a proton Q 2 s (x = 10 −2 ) = 0.11-0.14 GeV 2 using (3.27), which is consistent with (yet slightly smaller than) estimates based from fits to F 2 DIS data [50]. Note that the saturation scale in large nuclei and at smaller x considerably exceeds that in a proton, yieldingqL ∼1 GeV 2 , where the use of perturbative techniques is commonly assumed to be legitimate.

Scaling properties of heavy-quarkonium suppression
Before comparing the model predictions to the other available data, we discuss in this section the expected scaling properties of ψ suppression in the present model.
Let us first mention that the nuclear dependence of quarkonium suppression is often parametrised as a power law, where α is assumed to be independent of A. The power law is empirical. It can be inferred in the Glauber picture of ψ absorption in the nucleus, S abs ≃ exp −cst · A 1/3 , and using the approximation A 1/3 ≃ log A, which is accurate to the 10% level for 5 ≤ A ≤ 200. However, the Glauber picture of nuclear absorption is expected to hold when the ψ energy E in the nucleus rest frame is small enough, see Section 3.5.1. The heuristic law (4.1) has no reason to be valid at high E where the compact color octet QQ pair crosses the nucleus and hadronizes far beyond. We checked that the J/ψ suppression expected in our model does not follow the parametrization (4.1). To illustrate this, the typical values of α are found to vary by up to 10% depending on whether J/ψ suppression in p-W collisions (in the E866 kinematics) is normalized either to p-p or p-Be collisions. Clearly the attenuation factor R pA should be preferred to the effective power α when discussing nuclear suppression. We thus focus on R pA rather than on α throughout our study, and now discuss its scaling properties.
In the energy loss model of Gavin and Milana [20], quarkonium suppression exhibits an approximate x F scaling, i.e., R pA is a function of x F but independent of √ s. Indeed, assuming that the shape of dσ ψ pp /dx F is independent of √ s, and considering the limit In the present approach, however, the approximate x F scaling of quarkonium suppression is broken for several reasons: 1. The transport coefficientq, and therefore the function F loss in (4.2), depends explicitly on x 2 at small x 2 < x 0 , see (3.14). As we shall see, this effect is particularly important at LHC energies; 2. As discussed in Section 3.2, the slope of the p-p production cross section does depend on √ s (see Tables 1 and 2); 3. Finally, the saturation (or nPDF) effects also scale with x 2 yet this effect is actually rather small.
In order to illustrate this, J/ψ suppression has been computed in p-W collisions as a function of x F at NA3 ( Note that at LHC, the variation of R J/ψ pA with x F is extremely fast at very small x F . This strong dependence comes from the small-x behavior of the transport coefficientq(x), Eq. (3.14), together with the fast variation of x = x 2 with x F at |x F | 10 −2 . As we shall see in Section 4.5 the variation of R J/ψ pA with y ∼ ln x F is naturally much smoother. In order to check experimentally whether J/ψ suppression scales with x F , it would be crucial to measure J/ψ production in p-A collisions at RHIC and LHC at large x F , say x F 0.1, which is out of reach with the present apparatus. Such measurements could in particular shed light on the x dependence of the transport coefficientq(x).

Predictions and comparison to J/ψ data
Onceq 0 is determined from the fitting procedure described in Section 4.1, the x F dependence of the J/ψ quenching factor R ψ pA can be predicted in any target nucleus and at any center-of-mass energy for which the absolute p-p differential cross section dσ ψ pp /dx F has been measured. In this section we systematically compare the model predictions with all available data.

E866, NA3, E537, NA60, HERA-B
Let us start with the comparison of J/ψ suppression expected in an iron target and the E866 data for R Fe/Be , i.e., taken at the same energy as the fitted ratio R W/Be . The excellent agreement reported in Fig. 10 fully supports the atomic mass dependence of the model. This is at variance with the calculations by Gavin and Milana [20] which overestimated the ratio R Fe/Be , at that time measured by E772. 26 It is therefore a hint that the Ldependence expected here, ∆E ∝ √ L (see (1.3)), is probably more appropriate than the ad hoc assumption of Ref. [20], ∆E ∝ L.
Data taken at lower √ s or smaller x F are also compared to the model. As can be seen in Fig. 11 the agreement with NA3 p-A and π − -A data is excellent, both in shape and magnitude, over a very wide range in x F . It is also remarkable that the model is able to reproduce the different magnitude of suppression in p-A and π − -A collisions reported by NA3 [12]. This difference cannot be understood within nuclear absorption models, where nuclear suppression is a purely final state effect, thus independent of the projectile type. It cannot either be explained by nPDF effects, unless the nPDF to proton PDF ratios for valence quarks and for gluons, probed respectively in π − -A and p-A collisions, prove completely different. 27 In our picture, the smaller J/ψ suppression in π − -A collisions naturally arises from the flatter differential cross section, n πp = 1.4 vs. n pp = 4.3 at √ s = 19.4 GeV, see Table 1. Although no prediction of the exponent n is made in our model, it is clear that this feature can be explained from the larger slope, at large x, of the gluon PDF in a proton, xG(x) ∼ (1 − x) 3 [54], when compared to that of a valence antiquark PDF in a pion, xq(x) ∼ (1 − x) [55]. In this respect, let us mention that our assumption of an incoming gluon in quarkonium hadroproduction (see the Introduction) does not hold for NA3 pion-nucleus collisions, where subprocesses with an incoming valence antiquark dominate. In spite of this, a very good agreement between the model and the NA3 π − -A data is found, suggesting a mild dependence of the energy loss on the incoming parton type. A similar remark applies to the case of the π − -A E537 data discussed below. The E537 experiment also reported on measurements of J/ψ production in π − induced collisions on various nuclear targets (Be, Cu, W) at √ s = 15.3 GeV [56]. Our results 28 are found in reasonable agreement with the measured ratio R W/Be (Fig. 12, left); the slight underestimation of the suppression by the model might be attributed to nuclear absorption. Indeed, at this energy x crit F ≃ 0.7, and all E537 data lie in the x F ≤ x crit F domain. This might also explain the (more pronounced) difference between the observed and predicted magnitudes of the ratio R W/Cu (Fig. 12, right).
In Figs. 13 and 14 we compare our predictions with NA60 [43] 29 and HERA-B [32] p-A measurements. Although the center-of-mass energy is larger than those of NA3 and E537, the typical J/ψ energy range covered by NA60 and HERA-B is actually lower because of the smaller x F values probed by these experiments. As a consequence, J/ψ suppression can be affected more strongly by nuclear absorption effects, as can be inferred by the position of the x crit F arrows in Figs. 13 and 14, below which hadronization typically takes places inside the nuclear target. 30 26 We do not show the agreement between our model predictions and J/ψ E772 data [52] since those measurements were superseded by E866 [9]. 27 On top of this, nPDF effects in the NA3 kinematical domain, x2 ∼ 0.1-0.2, are known to be small for both valence quarks and gluons, see for instance the discussion in [53]. 28 Lacking π − -p data at E537 energy, we choose the exponent n = 1.4 (see Table 1). 29 Lacking p-p data at NA60 energies, we choose the exponent n = 4.3 (see Table 1). 30 In the left panel of Fig. 13, the arrow is not visible as x crit NA3 √s = 22.9 GeV (π beam) Figure 11. NA3 J/ψ suppression data [12] in p-A and π − -A collisions compared to the energy loss model. Nevertheless, the model predictions prove in very good agreement with data. In particular, the enhancement of J/ψ production observed by HERA-B at very negative x F , x F −0.2 (see Fig. 14) is well reproduced by the model. (The origin of R pA > 1 can be simply understood from the positive slope of dσ ψ pp /dx F in the target fragmentation region, x F < 0, see the HERA-B data in Fig. 6.) There is however room for J/ψ absorption with a cross section of a few millibarns, as suggested by the slight overprediction of R pA at NA60 precisely in the region where J/ψ absorption can no longer be neglected.

PHENIX
The predictions in d-Au collisions at RHIC, √ s = 200 GeV, are shown in Fig. 15 in comparison with PHENIX data [13], with (dashed line) and without (solid line) saturation effects. The energy loss model is able to reproduce nicely the J/ψ suppression at all rapidities. Note that in several phenomenological analyses, the suppression observed in the most forward rapidity bins has often been attributed to gluon saturation effects or to strong small-x shadowing in the nuclear PDF (see e.g. [47,57]). Here, the sole energy loss effects might be responsible for the observed suppression, although saturation might play a role as well. As a matter of fact, the agreement is better when saturation is included.
Remarkably, an excellent agreement is also observed in some negative y bins, for which nuclear absorption might also play a role (at least for y < y crit = −1.1). We shall discuss further these data in Section 4.7 when comparing to the predictions including nPDF effects. The J/ψ suppression has also been measured by PHENIX for various d-Au centrality classes [13] and more recently as a function of its transverse momentum [58]. Discussing these data would go beyond the scope of the present article and is left for future work.

Predictions and comparison to Υ data
The above comparison between J/ψ suppression data and our model predictions supports both the medium length and energy dependence of the model. The mass dependence of heavy-quarkonium suppression can be studied by investigating the suppression of Υ states in p-A collisions. Unfortunately the data are rather scarce; to our knowledge, the measurements have only been performed by E772 at Fermilab [59] and PHENIX and STAR [36,60] at RHIC.
The E772 data are shown in Fig. 16 for various nuclear targets (Ca, Fe, W) and in comparison to the model predictions. A rather good agreement between data and theory is found for x F > x crit F , although smaller experimental uncertainties would be necessary to further check the M dependence of the model. At small x F < x crit F , the measured nuclear production ratio R Υ pA lies much below our predictions, probably too much to be accommodated by Υ nuclear absorption. However, let us mention that the E772 measurements at low x F might be affected by uncorrected acceptance effects due to the correlation between x F and p ⊥ (see the discussion in [9]); it is therefore difficult to draw any conclusion from the significant disagreement observed at negative x F .  Figure 16. E772 Υ suppression data [59] in p-A collisions compared to the energy loss model.
These Υ data nevertheless allow the mass dependence of the energy loss to be constrained. In their paper [20], Gavin and Milana assumed that the mean energy loss scales as ∆E ∝ M −n , and considered explicitely the cases n = 2 ("power suppressed") and n = 0. From the comparison of their calculations with E772 data, these authors concluded that neither of these two choices were satisfactory: assuming ∆E ∝ M −2 led to too little Υ attenuation while a too strong suppression was predicted with the hypothesis ∆E ∝ M 0 . It is therefore interesting to note that the scaling ∆E ∝ M −1 predicted in [15] and used here (see (1.3)) supports this empirical observation, as the agreement in Fig. 16 indicates.
The predictions at RHIC are shown in Fig. 17. As expected, the suppression is less   [36,60], see Fig. 17. Hopefully more precise data will soon allow for clarifying the strength of Υ suppression in cold nuclear matter.

LHC predictions
We discuss in this section the J/ψ and Υ suppression expected in p-Pb collisions at the LHC ( √ s = 5 TeV). In Fig. 18 we show the R pPb ratios for both states as a function of the rapidity in the center-of-mass frame. 32 Interestingly, J/ψ production is significantly suppressed at large positive rapidity, e.g. R J/ψ pPb ≃ 0.7-0.8 at y = 1 and down to R pPb 0.5 at y 4. Because of the high center-of-mass energy of the collision at the LHC, saturation effects in the J/ψ channel are significant: in addition to energy loss, the suppression due to saturation ranges from S J/ψ A ≃ 0.9 in the most negative rapidity bins down to S J/ψ A ≃ 0.65 at y = 5. In the target fragmentation region (y < 0), the suppression is moderate (∼ 10-20%) while a possible J/ψ enhancement might be seen at very backward rapidities, y −5. Predictions using EPS09 [61] and DSSZ [62] nPDF sets are also discussed in Section 4.7.
In the Υ channel, the suppression is more moderate because of the mass dependence of energy loss, 33 ∆E ∝ M ⊥ −1 , e.g. R Υ pPb ≃ 0.85 at y = 3. At the LHC the saturation effects in the Υ channel prove quite small, although more pronounced than at RHIC. As can be seen from the arrow in Fig. 18, J/ψ and Υ hadronization should take place outside 31 The PHENIX and STAR data correspond to the suppression of Υ(1S)+Υ(2S)+Υ(3S) states. 32 Note that in p-Pb collisions at the LHC, the laboratory frame is shifted by ∆y ≃ 0.47 with respect to the center-of-mass frame. 33 Another effect, yet rather marginal, comes from the flatter x F distributions in p-p collisions (see Table 2 in Section 3.2). the nuclear medium for y > y crit ≃ −5; nuclear absorption should thus play little or no role at the LHC. These predictions can be compared to the future measurements during the p-Pb run scheduled at the LHC in January 2013. In order to test the model, the nuclear dependence of ψ production should ideally be measured for various rapidity bins and on a rather broad range, which hopefully should be possible with the ALICE or LHCb experiments. 34

E906 predictions
The E906 "SeaQuest" collaboration [63] aims at measuring Drell-Yan production in p-p and p-A collisions at E p = 120 GeV ( √ s = 15 GeV) at Fermilab. Although the first goal of this experiment is to study the sea quark asymmetry in the nucleon, it will also be able to measure the nuclear dependence of both Drell-Yan and J/ψ production on various nuclear targets and on a wide range in x F . In this section we present our model predictions on J/ψ suppression in p-A collisions at E906 energy, to be compared to the measurements that might already be available in 2013.
In Fig. 19 we plot the predictions in p-Fe (left) and p-W (right) collisions. 35 The suppression is very pronounced especially at large x F , for which however the J/ψ production cross section should be extremely small. 34 We recall that the calculations are made using p ⊥ = 1 GeV in the calculation of the transverse mass. Therefore our predictions on J/ψ suppression should be adapted for the CMS acceptance which requires a transverse momentum cut, p ⊥ 6 GeV, in the J/ψ channel. 35 Lacking p-p data at this energy, we choose the exponent n = 4 to be consistent with the systematics discussed in Section 3.2.

Comparing predictions using saturation vs. nPDF
For completeness, we compare in this section the former results on J/ψ suppression at RHIC and LHC obtained in the "energy loss + saturation" model with the predictions using the EPS09 [61] and DSSZ [62] nPDF leading-order sets instead of saturation. Unlike saturation effects, nPDF corrections should be valid (and possibly non-negligible) even at not too small values of x 2 , and in particular at E866 energy. Therefore, the transport coefficientq 0 using each of the two nPDF sets has been consistently refitted to E866 data.
The corresponding values, used for the RHIC and LHC predictions with these two sets, are indicated in Table 4. The comparison is shown in Fig. 20 at RHIC. The predictions using saturation, EPS09 and DSSZ somehow differ in the rapidity dependence of R J/ψ dAu . The DSSZ modifications are rather small, leading to a suppression similar to the one assuming energy loss effects only. On the contrary, the EPS09 set exhibits larger modifications to the gluon nPDF (and in particular a slightly faster variation in this x domain) increasing the slope of R J/ψ dAu versus y. At mid-and forward rapidity, the various predictions are similar; in particular all of them reproduce nicely the data at y = 1-2 with a slightly better description with saturation or using the EPS09 set as compared to DSSZ (yet this is not statistically significant). On the other hand, the predictions are different at backward rapidity. The best agreement is obtained assuming saturation effects (in addition to parton energy loss) or DSSZ nPDF instead of EPS09. This observation obviously depends on the present energy loss model, preventing us from drawing a firmer conclusion.
In order to analyze a bit more quantitatively these results, the values of χ 2 /ndf obtained for the E866 and PHENIX data sets are quoted in Table 4. As can be seen, the E866 data do not allow one to differentiate the various energy loss predictions including (or not) EPS09/DSSZ nPDF corrections. On the contrary, the agreement at RHIC is considerably better when energy loss is supplemented by saturation effects (χ 2 /ndf = 0.3) rather than  by nPDF (χ 2 EPS09 /ndf = 2.7, χ 2 DSSZ /ndf = 1.1), as mentioned above. For completeness we also quote the values of χ 2 /ndf assuming no energy loss but only saturation 36 or nPDF corrections. As can be seen from Table 4, saturation without energy loss still gives a fair description of PHENIX data, χ 2 /ndf = 2.7 (as well as EPS09 and DSSZ to a lesser extent), although saturation and nPDF effects alone would totally fail to reproduce E866 data. 36 The value ofq0 = 0.075 GeV 2 /fm quoted in Table 4 is here to determine the saturation scale. q 0 (GeV 2 /fm)  Table 4.q 0 and χ 2 /ndf of E866 and PHENIX data, with (upper rows) and without (lower) energy loss effects, for various assumptions regarding the nuclear modifications of gluon distributions in nuclei.
The predictions at the LHC are shown in Fig. 21 (left). The expected suppression with saturation effects or using the EPS09 set prove remarkably similar. On the contrary, the nPDF corrections given by DSSZ are tiny in the forward rapidity bins, despite the small values of x probed in the Pb nucleus. Although the absolute magnitude of the J/ψ suppression differs depending on the assumption regarding saturation/nPDF effects, the rapidity dependence (especially at y > 0) is mostly governed by energy loss effects. This could be tested experimentally. In the present model, energy loss effects are moderate at mid-rapidity which corresponds to the maximum of the p-p production cross section. As a consequence, the expected suppression at y = 0, R pPb (y = 0), is more sensitive to saturation/nPDF effects. Moreover, since the rapidity dependence is essentially due to energy loss effects, the double ratio R pPb (y)/R pPb (y = 0) is rather independent of the strength of saturation/nPDF effects. This is illustrated in Fig. 21 (right) where R pPb (y)/R pPb (y = 0) is plotted. As can be seen this double ratio proves remarkably similar whether or not energy loss is supplemented with nPDF or saturation effects.

Discussion
The agreement between our model and the p-A data for quarkonium nuclear suppression is quite remarkable. With a single free parameterq 0 , both the slope and normalization of R pA (or R pA /R pB and also R πA ) are accurately described, for various collision energies, various target nuclei and different masses (J/ψ, Υ), and over a broad range in x F (or rapidity). We also stressed that the effect of saturation alone fails in describing J/ψ nuclear suppression at different collision energies. This strongly supports parton energy loss as a dominant effect in p-A quarkonium nuclear suppression, the main conclusion of our study. The successful description of the data is mostly due to the (medium-induced) energy loss scaling as ∆E ∝ E, where E is the energy of the QQ pair in the nucleus rest frame. This behaviour arises when the partonic subprocess looks like small angle scattering of an asymptotic charge, and thus holds within our assumption of a long-lived, color-octet QQ pair. Our results support the parametric (E, M and L) dependence of the induced radiation spectrum (1.2), which is derived from first principles in Section 2.
These results also give some hint on the mechanism for low p ⊥ (p ⊥ M ) heavyquarkonium hadroproduction. We argued in the Introduction that at large x F , the octet QQ pair should be long-lived in any quarkonium production model, including the Color Singlet Model (CSM). The agreement of the energy loss model with the large x F suppression data thus cannot distinguish between production models. But we found that the agreement extends to small values of x F (see in particular the comparison to RHIC data at y ∼ 1-2, corresponding to x F ∼ 0.04-0.1, in Fig. 15), where assuming a long-lived coloroctet QQ pair becomes inaccurate in the CSM. The CSM mechanism thus seems somewhat disfavoured by our results, at least as a dominant contribution to inclusive (i.e., low p ⊥ and low x F ) J/ψ production. The future measurements in p-Pb collisions at the LHC will probe small values of x F (|x F | < 0.1), yet in a rather large rapidity interval (|y| < 5), and might thus further clarify the underlying dynamics of heavy-quarkonium production. In fact our energy loss explanation of J/ψ suppression is consistent with any J/ψ production model where t hard ≪ t octet , leaving room for gluon radiation with t hard ≪ t f ≪ t octet , see (2.18). It was argued in Ref. [30] that a qualitative analysis of the quarkonium production data suggests a mechanism for quarkonium hadroproduction, named "Comover Enhancement Scenario", where color neutralization is realized at the time t octet by a semi-hard scattering between the QQ pair and the comoving radiation field of the incoming parton. It is intriguing that the condition on t octet inferred from global production features, t hard ≪ t octet ≪ t ψ [30], is consistent with the condition (2.18) necessary to explain nuclear suppression from radiative parton energy loss.
The induced radiation spectrum was derived assuming a given hierarchy between the nuclear size L, gluon formation time t f and quarkonium hadronization time t ψ . Thus, as we emphasized several times, the model should in principle be valid only when the quarkonium state hadronizes outside the nucleus, i.e., when E is large enough or x F > x crit F . It is quite striking that the extrapolation of the model to the region x F < x crit F is either consistent with the data (within error bars, see e.g. the NA3 data in Fig 11, HERA-B data in Fig. 14 and PHENIX data in Fig. 15), or systematically underestimates quarkonium nuclear suppression (NA60 data in Fig. 13). This suggests parton energy loss to play a role in a broader domain than expected, the possible additional suppression required at x F < x crit F being due to nuclear absorption of the fully formed quarkonium state. In our study we also assumed quarkonium production in proton-nucleus collisions to arise from the splitting of an incoming gluon into an octet QQ pair. This assumption becomes inaccurate at very large x F , where quark-induced processes (such as qq → QQ) come into play. Although we expect the parametric dependence of the associated radiation spectrum to be unchanged, the overall color factor might be changed in this case. A possibly smaller effective color factor at very large x F might explain the milder J/ψ suppression observed by E866 at x F ∼ > 0.8 (see Fig. 8) than predicted in our model. However, as already mentioned in Section 4.3.1, the very good agreement between the model and the NA3 pion-nucleus data, for which the qq annihilation channel is dominant at all x F , suggests a relatively weak dependence of the energy loss on the incoming parton type.
We might envisage refinements of the parton energy loss model presented here, such as including quarkonium absorption at x F < x crit F and quark-induced processes, in order to extend the domain of validity of our approach. However, we find it more important to first confirm the dominant role of parton energy loss in p-A collisions, where gluon-induced processes dominate and our model assumptions apply.
First, the energy loss model can be tested in forthcoming p-Pb collisions at the LHC, for which our predictions for the y-dependence of J/ψ and Υ suppression are shown in Fig. 18, and in p-A collisions in the E906 fixed-target experiment at Fermilab (Fig. 19). The model should as well be confronted to the existing RHIC d-Au data on the p ⊥ and centrality dependence of J/ψ suppression, measured at various rapidities [13,58]. This requires generalizing (3.6) to double differential (in x F and p ⊥ ) cross sections, and will be the subject of a future work. It will be interesting to check whether the L-dependence of the energy loss predicted in (1.3) is consistent with the centrality dependence of the RHIC d-Au data. Our study might also help interpreting quantitatively quarkonium measurements performed in heavy-ion collisions at RHIC [64,65] and LHC [66][67][68]. Indeed, parton energy loss through the incoming cold nuclei is expected to combine with hot effects (such as Debye screening or final state energy loss in a QGP). The evaluation of J/ψ suppression in A-A collisions expected from cold parton energy loss alone will be presented in a future study. Finally, as discussed in the Introduction other processes than quarkonium production should be sensitive to a parametrically similar parton energy loss, such as open charm and light hadron production in p-A collisions. This work is also in progress.

A x dependence ofq
In Ref. [38], the transport coefficientq was related to the gluon distribution G(x) in a target nucleon, see the expression (3.13). The value of x to be used in xG(x) in (3.13) can be estimated by considering the kinematics of the rescattering process.
Following Ref. [38], let us consider the specific case of DIS, where an energetic light quark of momentum p is produced and then rescatters with transfer q on a target nucleon of momentum P . Working in the target nucleus rest frame, we have P = (m N , 0), with m N the nucleon mass. Choosing light-cone coordinates p ± = p 0 ± p z and p = (p + , p − , 0 ⊥ ) along the negative z-direction, the condition for the final quark to be on-shell reads (p + q) 2 = (p + + q + )(p − + q − ) − q 2 where we neglected q − as compared to p − = 2E.

parton produced inside the target
When the hard production time t hard ∼ E/Q 2 ≪ L, or equivalently when the Bjorken variable x B ≡ Q 2 /(2m N E) ≫ x 0 ≡ 1/(2m N L), the parton p is effectively produced incoherently inside the nucleus. This is the situation considered in Ref. [38], which we now briefly review. If the rescattering occurs at a distance z from the production point, from the uncertainty principle we have |p + | ∼ 1/z just before the scattering. From the constraint z ≤ L, we obtain |p + | ∼ > 1/L. For p − large enough Eq. (A.1) gives q + ≃ |p + |. The momentum fraction of the rescattering gluon thus satisfies [38] x ≡ q parton produced far before the target When t hard ≫ L ⇔ x B ≪ x 0 , the virtual photon splits into a light quark-antiquark pair far before the nucleus, and the DIS process is coherent over the whole nucleus. In this case, the quark virtuality |p 2 | = |p + p − | ∼ Q 2 , and |p + | ∼ Q 2 /p − is not bounded by 1/L any longer. From Eq. (A.1) we obtain (using q 2 ⊥ ≪ Q 2 ) For a generic hard process (for instance in a p-A collision) of coherence length t hard ∼ E/M 2 ∼ 1/(2m N x 2 ), the above DIS example supports the following estimate for the value of x to be used in Eq. (3.13), x = x 0 Θ(x 2 > x 0 ) + x 2 Θ(x 2 < x 0 ) = min(x 0 , x 2 ) ; x 0 ≡ 1 2m N L , (A.4) thus specifying the x 2 -dependence of the transport coefficientq.