Centrality and pT dependence of J/psi suppression in proton-nucleus collisions from parton energy loss

The effects of parton energy loss and pT-broadening in cold nuclear matter on the pT and centrality dependence, at various rapidities, of J/psi suppression in p-A collisions are investigated. Calculations are systematically compared to E866 and PHENIX measurements. The very good agreement between the data and the theoretical expectations further supports pT-broadening and the associated medium-induced parton energy loss as dominant effects in J/psi suppression in high-energy p-A collisions. Predictions for J/psi (and Upsilon) suppression in p-Pb collisions at the LHC are given.


Introduction
A wide range of phenomena observed in heavy-ion collisions at RHIC and LHC, both in hard and soft processes, suggests that a new, strongly interacting state-of-matter has been created, the quark-gluon plasma (QGP). One of the most striking observables is the suppression of high-p ⊥ particles ("jet-quenching") in A-A compared to a simple scaling of p-p collisions [1][2][3][4], which together with dijet momentum imbalance [5,6] finds a natural explanation in parton energy loss in the QGP [7]. Another observable in heavy-ion collisions which received much attention is heavy-quarkonium suppression. Such a suppression is expected from the Debye screening of the in-medium heavy-quark potential, and was thus originally proposed as a potential signal of QGP formation (and a direct probe of the plasma temperature) [8]. However it was later realized that several other effects can modify quarkonium yields in A-A collisions (see [9] for a review), some of those effects playing a role even in absence of any hot medium.
In order to quantify the properties of the QGP created in heavy-ion collisions -a main goal of the forthcoming measurements at RHIC-II and LHC -a solid understanding of the nuclear modification of particle spectra in cold nuclear matter is thus required. A baseline for the study is provided by p-A (or d-A) collisions where a significant suppression is already reported for some particle species. In particular, light hadron [10,11] and J/ψ [12,13] production in p-A collisions at forward rapidity is significantly lower than expected from the naive scaling of p-p spectra.
At the same time it is astonishing that no consensus on the cold nuclear matter effects responsible for J/ψ suppression has been achieved yet. In addition to the modifications of nuclear parton distribution functions, various mechanisms have been proposed to explain J/ψ suppression in p-A collisions. In the nucleus rest frame, a high-energy J/ψ is formed long after the nucleus thus what actually propagates through the nucleus is the parent cc pair. Some approaches attribute J/ψ suppression to an effective absorption cross section σ abs of the cc pair (see Refs. [14] and [15] for recent works). Other models attribute J/ψ suppression to the increase of the cc pair invariant mass by the multiple soft rescatterings through the nucleus, leading to a reduction of the overlap with the J/ψ wave function [16] (see also [17]). In the approach used by [18] and that of our group [19], the dominant role is played by parton radiative energy loss, i.e., gluon radiation induced by multiple scattering of fast partons (or color octet cc pair) travelling through the nucleus.
The present study, together with our previous works [19][20][21], supports parton energy loss induced by p ⊥ -broadening as the main effect in J/ψ suppression. At this point we should stress that the physical content of our approach and that of [18], both based on parton energy loss, are actually quite different. In [18], it is assumed that the interference between gluon emissions off the incoming and outgoing parton participating to the hard production process can be neglected. Under this assumption the induced energy loss of an incoming gluon in J/ψ production is parametrically the same (up to color factors) as that of the incoming quark in the case of Drell-Yan pair production [22].
It was however argued in [21] that when the outgoing parton (or compact cc pair in the case of J/ψ production) is produced at small angle (i.e., large energy E at limited p ⊥ ) in the target nucleus rest frame, the medium-induced gluon spectrum is dominated by the interference between initial and final state radiation of gluons with large formation times. Such radiation is expected in low-p ⊥ J/ψ production in p-A collisions, where an incoming gluon is scattered at small angle into a compact color octet cc pair. The associated energy loss is proportional to the J/ψ energy E and the role of this effect is thus expected to increase with increasing rapidity. This observation is at the basis of the energy loss scenario used in [19] to describe J/ψ suppression as a function of y (or x F ) in a wide collision energy range for minimum bias p-A and d-A collisions. In the present study we generalize this approach to address the p ⊥ and centrality dependence of nuclear suppression at fixed-target (E866), RHIC and LHC energies.
The paper is organized as follows. In Section 2 we generalize the calculation scheme of [19] to address the p ⊥ and centrality dependence of J/ψ nuclear modification factors. Section 3 is devoted to the comparison of the model to E866 and latest PHENIX data and predictions in p-Pb collisions at the LHC are given in Section 4. Results are summarized in Section 5.

Shift in E and p ⊥
Heavy-quarkonium 1 nuclear suppression in minimum bias p-A collisions as compared to p-p collisions can be represented in terms of the ratio where y and p ⊥ are the quarkonium rapidity and transverse momentum in the c.m. frame of an elementary p-N collision (of energy √ s). By convention forward (positive) rapidities correspond to the proton fragmentation region.
In the present study we generalize the model of Ref. [19] by expressing the quarkonium double differential cross section in p-A collisions in terms of that in p-p collisions, where a shift in the quarkonium energy E (defined in the nucleus rest frame) as in Ref. [19] but also a shift in p ⊥ account, respectively, for the energy loss ε and transverse momentum broadening ∆ p ⊥ of the octet QQ pair propagating through the nucleus. The "double shift" in E and p ⊥ relates the p-A and p-p double differential cross sections dσ/dEd 2 p ⊥ as 2 The quantity P(ε, E) is the energy loss probability distribution or quenching weight (see Section 2.2) associated to the medium-induced radiation spectrum in a target nucleus A as compared to a lighter (e.g., proton) target nucleus. The integral over ε is bounded by is the projectile proton energy in the nucleus rest frame. (We work in the limit √ s ≫ m p , with m p the proton mass.) For the time being we consider minimum bias p-A collisions, hence the normalization factor 1/A in the l.h.s. of (2.2), but the model will be generalized in Section 2.5 to p-A and d-A collisions in a given centrality class. We will assume that ∆ p ⊥ is uniformly distributed in the azimuthal angle ϕ, and has a modulus defined by (2.10). The relation between the p-A and p-p differential cross sections in y and p ⊥ , can be simply obtained from (2.2) by using

4)
1 denoted by "ψ" in the rest of the paper. 2 The dependence of ∆ p ⊥ on the azimuthal angle ϕ, ∆ p ⊥ ≡ ∆ p ⊥ (ϕ), will be implicit in the following, and we will use the notations: where M ⊥ = (M 2 + p 2 ⊥ ) 1 2 , with M the mass of the QQ pair. In Eq. (2.3), the variable E is given by E = E (y, p ⊥ ) = E p e y M ⊥ / √ s, and the p-p cross section is evaluated at For the double differential p-p cross section dσ pp /dyd 2 p ⊥ , we will use a parametrization consistent with the p-p data, see Section 2.3. In our model, the nuclear modification factor (2.1) obtained from (2.3) is thus fully determined by the quenching weight P(ε, E) and the transverse momentum broadening ∆p ⊥ .

Quenching weight and ∆p ⊥
In Ref. [19] the appropriate medium-induced gluon radiation spectrum dI/dω in quarkonium production in p-A collisions (as compared to p-B collisions) was derived, as well as the associated quenching weight P(ε, E), where ℓ ⊥A represents the accumulated transverse momentum transfer due to soft rescatterings in the target nucleus A, and Λ B = max(Λ QCD , ℓ ⊥B ). The latter dependence on Λ QCD arises from vetoing induced gluon radiation with k ⊥ < Λ QCD [19]. Note that P(ε, E) is determined analytically in terms of the dilogarithm function Li 2 (x), (2.7) The semi-hard transfer ℓ ⊥A is given by where L A is the effective path length 3 across the target nucleus A andq A the transport coefficient [19] The parameterq 0 = 0.075 GeV 2 /fm was extracted in Ref. [19] from a fit to the E866 data for R J/ψ W/Be . We will use this value in the present study, which thus contains no free parameter. 4 3 For the present study we will use LBe = 3.24 fm, LFe = 6.62 fm, LW = 9.35 fm, LAu = 10.23 fm and L Pb = 10.11 fm for minimum bias collisions, as well as Lp = 1.5 fm [19]. The dependence of LAu (resp. L Pb ) on the centrality class in d-Au collisions at RHIC (resp. p-Pb collisions at LHC) is described in Section 2.5, see Table 2, and in Appendix B. 4 Let us also note that in the fitting procedure of Ref. [19], small values ofq0 and thus ofq L were explored, for which the dependence of (2.6) on Λ QCD is relevant. However, it turns out that the extracted q0 = 0.075 GeV 2 /fm is large enough to satisfy √q L > 0.25 GeV for all values of LA (including Lp) and x2 considered in the present study, as can be easily checked from (2.9) using Lp = 1.5 fm. The present study is thus independent of the value of Λ QCD provided Λ QCD ≤ 0.25 GeV.
Finally, the transverse momentum broadening in p-A with respect to p-B collisions is simply given by (2.10) In our study we neglect the fluctuations of the broadening ℓ 2 ⊥ around the average valuê qL. We have checked that assuming a distribution in ℓ ⊥ of the Gaussian type, P (ℓ ⊥ ) ∝ exp (−ℓ 2 ⊥ /qL), only slightly modifies (by at most 10%) the predictions for R ψ pA (y, p ⊥ ) presented in Sections 3 and 4, without changing the overall shape of R pA .

Parametrization of the p-p cross section
Similarly to Ref. [19], we use for the double differential p-p cross section dσ pp /dyd 2 p ⊥ a simple parametrization consistent with the available p-p data, rather than relying on some model-dependent quarkonium production mechanism in hadronic collisions.
The double differential cross section of prompt J/ψ and Υ production can be conveniently parametrized as At the LHC ( √ s = 7 TeV) the values of the free parameters N , p 0 , m and n are obtained from a global fit of ALICE [23], ATLAS [24] and LHCb [25] data on prompt J/ψ production 5 and from a fit of LHCb data [26] on Υ production. They are summarized in Table 1 together with the corresponding At RHIC the amount and precision of data are not sufficient to fix precisely the fit parameters. The fit to the double differential data measured by PHENIX [27] is therefore performed by fixing the value of n = 8.3 obtained from the fit of the single differential cross section dσ J/ψ pp dy, see Ref. [19]. Finally, a fit to the J/ψ E789 ( √ s = 38.7 GeV) data [28] is also performed with the aim to compare the model predictions to the E866 data [12] at the same center-of-mass energy. 7 Like at RHIC, the number of fit parameters is reduced by fixing the value of n = 4.5 obtained from fitting dσ J/ψ pp dx F data at this energy [19]. B µµ dσ/dp T (pb/GeV) LHCb √s = 7 TeV 2 ≤ |y| ≤ 2.5 ϒ Figure 1. Comparison between heavy-quarkonium (J/ψ, Υ) production data in p-p collisions and the parametrization (2.11) (solid red line). Data are taken from ATLAS [24], LHCb [25,26], and PHENIX [27].
The agreement between some of the prompt J/ψ and Υ measurements at RHIC and LHC and the parametrization (2.11) is shown in Fig. 1. Finally, let us mention that the functional form of the parametrization (2.11) is fully consistent with that for the single differential rate in x F used in [19], as we briefly show in Appendix A.

A useful approximation for R ψ
pA (y, p ⊥ ) Using (2.3) and (2.11) the attenuation factor (2.1) reads Since P(ε, E) is peaked at small values of ε, we neglect ε in the latter ϕ integral. In this approximation, the ϕ and ε integrals factorize, The factor R broad pA (y, p ⊥ ) describes nuclear modification due to transverse momentum broadening only, as can be seen by setting P(ε, E) = δ(ε) in (2.12). The factor R loss pA (y, p ⊥ ) describes the effect of energy loss only, obtained by setting ∆p ⊥ = 0 in (2.12). In the following we will use the factorized expression (2.13), which turns out to be a very accurate approximation to (2.1) in all the practical applications of Sections 3 and 4.
Finally, let us mention that the p ⊥ -inclusive nuclear suppression factor R ψ pA (y) studied in Ref. [19] can be simply recovered from (2.12), along the same lines as in Appendix A where the parametrization of dσ pp /dy is obtained from dσ pp /dyd 2 p ⊥ . Integrating both the numerator and denominator of (2.1) (and thus of (2.12)) over p ⊥ , the function ν(y, p ⊥ ) can be replaced by its value at a typicalp ⊥ determined by the width of µ(p ⊥ ) (see (A.6)), and the p ⊥ -integral of µ(p ⊥ ) cancels between the numerator and denominator. As a result, the p ⊥ -inclusive suppression factor reads R ψ pA (y) ≃ R loss pA (y,p ⊥ ), which corresponds exactly to the quantity studied in Ref. [19].

Centrality dependence
In the preceding sections we addressed the case of minimum bias p-A collisions. The model is generalized to the case of p-A collisions at a given centrality (or at a given impact parameter b) by using the effective length L A corresponding to that centrality. Since L A fully determines the essential quantities ℓ ⊥A and ∆p ⊥ of the model (see (2.8), (2.9) and (2.10)), this is the only modification required. However, the impact parameter (as well as the number of participants) is not a direct experimental observable, and the consistent way for making a theoretical estimate of L A is to follow the experimental procedure as closely as possible.
Both at RHIC [29] and the LHC [30] a centrality selection is done by triggering on event multiplicity in forward detectors. This multiplicity is strongly correlated with the number of participating nucleons from the target nucleus. Thus, multiplicity cuts impose a restriction on the number of participants in a given event. The values of the forward event multiplicities which separate centrality classes in the experiment are chosen to attribute a certain fraction of the total inelastic cross section to each centrality class. The common choice is making 20% slices from the most central (largest multiplicity) to the peripheral (lowest multiplicity) events. The exception is the most peripheral class which is taken at 60-88% of the total inelastic cross section at RHIC (class D) and 60-100% at the LHC (class 4).
To compute the average path length L A for each centrality class, we employ the following procedure.
First, we define centrality classes in terms of the number of participants N p . Within a Glauber description (see Appendix B.1), we define a centrality class by the threshold values N min p and N max p of N p which saturate approximately the same fraction of the total interaction probability P (N p ≥ 1) = 1 as the fraction of the total inelastic cross section attributed to the centrality classes in the experimental selection procedure. We note however that centrality in the experiment is defined in terms of event multiplicity rather than N p . In a given multiplicity class one may have events where N p is slightly above or below the thresholds defined according to a sharp cut on the fraction of the total interaction probability as a function of N p . In order to account for this possibility we widen the N p interval attributed to each class and use the new threshold values in further estimates. Those values are given in Table 2. As a consistency check of our class selection method, we calculate the average number of binary collisions N c for the different centrality classes. These numbers coincide with the results of the Glauber Monte-Carlo supplemented with the p-p RHIC data used by the PHENIX collaboration (see Table 1 in Ref. [13]).
Second, within each of the centrality classes defined in such a way, we determine, also in the Glauber model, the average number of target nucleons participating in the rescattering of the fast color octet QQ pair, see Appendix B.2. The average path length L A for a given centrality class directly follows from this number. Table 2 displays the values of L Au (RHIC) and L Pb (LHC) to be used in our model, for each centrality class.

Comparison to E866 and RHIC data
The J/ψ nuclear production ratio is computed as a function of p ⊥ for various values of y or x F = (2M ⊥ / √ s) sinh y. The only parameter of the model, the transport coefficientq 0 , has been fixed toq 0 = 0.075 GeV 2 /fm in [19] from the x F dependence of J/ψ suppression measured in p-W collisions by E866 [12]. In the numerical calculations we use Λ QCD = 0.25 GeV (however see footnote 4), α s = 0.5, and M = 3 GeV (M = 9 GeV) for the mass of the cc (bb) pair.

E866
The E866 collaboration has measured the J/ψ suppression in p-Fe and p-W collisions (with respect to p-Be) at √ s = 38.7 GeV as a function of the transverse momentum for three domains in x F [12]. At small x F , J/ψ production can be affected by nuclear absorption since the typical J/ψ hadronization time becomes comparable to (or less than) the size of the nuclear targets [19].  In Fig. 2 are shown as solid lines the Fe/Be (left) and W/Be (right) nuclear production ratios at intermediate (top) and large (bottom) x F . The ratio R pA (p ⊥ ) increases with p ⊥ in almost the whole p ⊥ range, with more pronounced suppression (at p ⊥ = 0 GeV, where R pA is the smallest) in W targets and at large x F . The p ⊥ dependence essentially arises from that of R broad pA (p ⊥ ), shown as dashed lines. The sole energy loss effect, R loss pA (y, p ⊥ ), proves rather flat in this p ⊥ -domain, but is essential to fix the magnitude of R pA , leading to a remarkable agreement between the data and the model predictions in the 0 ≤ p ⊥ ≤ 2 GeV range.
At large x F , the data overshoot the theoretical expectations above p ⊥ 3 GeV, for which the model predictions flatten out. Several reasons might explain this disagreement. First of all, the parametrization of the p-p cross section (using E789 data, see Section 2.3) has been performed at x F = 0 only; it could therefore well be that this fit, used to compute R pA , is no longer appropriate to describe the p-p production cross section at both large x F ≃ 0.4 and p ⊥ 3 GeV. Also note that the theoretical calculations have been performed at x F = x F while the data are averaged on a rather large x F -bin. When p ⊥ gets larger, the kinematical correlation between p ⊥ and x F would therefore tend to decrease the typical x F , for which the model would predict slightly larger R pA ratios. 8 For those reasons it is difficult to draw any firm conclusion from the comparison between the model and the E866 data at large x F and p ⊥ .
Putting aside the latter (p ⊥ , x F ) region, it is remarkable that the model reproduces quantitatively the p ⊥ dependence of J/ψ suppression in various nuclei and at different x F values. The fact that the same quantity,qL, determines the strength of medium-induced gluon radiation (and therefore energy loss) necessary to explain the x F dependence of R pA (see Ref. [19]) and the amount of momentum broadening required in the present study to reproduce the shape of R pA as a function of p ⊥ , strongly supports parton energy loss induced by momentum broadening as the dominant effect in quarkonium nuclear suppression.

RHIC
Let us now move to RHIC energy, where the p ⊥ dependence of J/ψ suppression in d-Au collisions has been reported recently by the PHENIX collaboration [13].
In Fig. 3 the model predictions are compared to the PHENIX data measured in minimum bias d-Au collisions, at backward (−2.2 ≤ y ≤ −1.2, left), central (|y| ≤ 0.35, middle) and forward (1.2 ≤ y ≤ 2.2, right) rapidities. 9 The model reproduces the trend seen in data, namely an increase with p ⊥ up to p ⊥ ≃ 4-5 GeV. Around those values, some nuclear enhancement, R pA (p ⊥ ) > 1, is visible at backward and central rapidities, but would need more precise data to be confirmed. At forward rapidity, the suppression due to energy loss is too strong to observe such an enhancement, both in the data and in the model. As for the results at E866 energy discussed earlier, the p ⊥ shape of R pA (p ⊥ ) is essentially driven by the effect of momentum broadening, Eq. (2.14), shown as dashed lines.
On top of minimum bias collisions, J/ψ suppression has also been measured in four centrality classes of d-Au collisions (0-20%, 20-40%, 40-60%, 60-88%). The data are shown in Fig. 4 for the three rapidity bins and four centrality classes and compared to the model. Details on the medium length L corresponding to the various centrality classes of d-Au collisions and used in the theoretical calculation can be found in Section 2.5 and Appendix B. As for minimum bias collisions, the model predictions prove in very good agreement with data. In particular, the centrality dependence is well reproduced by the model, with rather pronounced effects in the 0-20% most central collisions, and almost negligible effects, R pA ≃ 1 in the whole p ⊥ -range for the 60-88% most peripheral collisions.
Therefore it appears that the present model, based on parton energy loss and momentum broadening, offers a better agreement on the p ⊥ and centrality dependence of R pA than models based purely on modification of parton densities [15] and nuclear absorption [14] which tend to predict a flatter dependence of R pA (p ⊥ ).

Predictions for p-Pb collisions at the LHC
A run of p-Pb collisions at √ s = 5 TeV has taken place at LHC early 2013 with an integrated luminosity of L ≃ 30 nb −1 collected by ALICE, ATLAS and CMS experiments. This will allow for precise measurements of J/ψ and Υ production in p-A collisions at an unprecedented energy, on a wide range of rapidities (|y| < 5) and transverse momenta, and consequently clarify the role of cold nuclear matter effects at high energy. In this section, we therefore provide predictions for J/ψ and Υ nuclear production ratios R pA as a function of p ⊥ , for different rapidities (y = −3.7, 0, 2.8) 10 and centrality classes (labelled 1. . . 4 from central to peripheral collisions, see Section 2.5 for details).
In Fig. 5 we show the J/ψ nuclear production ratio for minimum bias p-Pb collisions, from backward (left) to central (middle) and forward (right) rapidities. At all rapidities a depletion of J/ψ production is expected at low p ⊥ , say p ⊥ 3 GeV. At larger p ⊥ a "Cronin peak" might only be visible at large rapidity, y = 2.8 and above, in minimum bias collisions. Quite generally the effects of both energy loss and momentum broadening are expected to become more pronounced at larger rapidities. As shown in [19], energy loss effects prove stronger at large positive rapidity because of the energy dependence of the average energy loss, ∆E ∝ E, associated to the medium-induced spectrum (2.6). LHC √s = 5 TeV y=2.8 Figure 5. Prediction for the J/ψ nuclear suppression factor R pA (p ⊥ ) in minimum bias p-Pb collisions at LHC, for backward, central and forward rapidities.
The model predictions are also provided in the four centrality classes of p-Pb collisions in Fig. 6. As expected the deviations of R pA from unity are largest in the most central collisions, while in the most peripheral p-Pb collisions (centrality class 4), R pA (p ⊥ ) ≃ 1 at all p ⊥ 2-3 GeV. The most spectacular effects can be seen in the central collisions (class 1) and at forward rapidity, where R pA ≃ 0.25 at p ⊥ = 0 GeV and R pA ≃ 1.3 at p ⊥ = 6 GeV. Predictions are performed as well in the Υ channel. The expected R Υ pA at the LHC is shown in Fig. 7 for the most central p-Pb collisions (centrality class 1). Because of the mass dependence of the energy loss, ∆E ∝ M −1 ⊥ , the suppression due to energy loss, R loss pA , is milder than for J/ψ. Moreover, although the amount of momentum broadening experienced by J/ψ and Υ states is expected to be similar (neglecting the x dependence of the transport coefficient), the height of the Cronin peak is much less pronounced for Υ than for J/ψ (compare e.g. Fig. 6 bottom and Fig. 7 right). 11 This is due to the flatter Υ p ⊥ -spectrum as compared to that of J/ψ production, see Fig. 1 and the values of the parameters p 0 and m in Table 1. As can be seen in Fig. 7, Υ suppression, R Υ pA < 1, is predicted in the range 0 ≤ p ⊥ ≤ 6 GeV at mid-rapidity. The suppression extends to larger p ⊥ than for J/ψ due to the larger value of the p 0 parameter in the p-p cross section. The suppression is maximal at p ⊥ = 0 GeV, where R Υ pA ≃ 0.85 (resp. 0.65) at y = 0 (resp. y = 2.8). 11 As a matter of fact, the slight enhancement arising from R broad pA is compensated by energy loss effects, R loss pA < 1, making R Υ pA smaller than one at all p ⊥ .

Conclusion
Following our earlier work [19], we studied the effects of parton p ⊥ -broadening and energy loss in cold nuclear matter on the p ⊥ dependence of J/ψ and Υ suppression in p-A collisions. We found that the momentum broadening is responsible for the fast variation of J/ψ suppression with p ⊥ , while medium-induced energy loss essentially affects the magnitude of R pA . Using the transport coefficientq 0 = 0.075 GeV 2 /fm fixed in [19], the model predictions prove in very good agreement with recent PHENIX data [13] in minimum bias and centrality-dependent d-Au collisions at √ s = 200 GeV. Our results are also successfully compared to earlier results from the E866 collaboration [12]. Finally, predictions for J/ψ and Υ suppression in p-Pb collisions (minimum bias and in four centrality classes) at the LHC ( √ s = 5 TeV) are provided.
The good description, within a consistent framework, of both the rapidity and transverse momentum dependence of R J/ψ pA from fixed-target experiments to RHIC is a hint that parton energy loss induced by momentum broadening might be the dominant effect responsible for J/ψ suppression in p-A collisions.

A Fits of p-p single and double differential rates
Here we compare the parametrization (2.11) for the double differential p-p cross section used in the present study and that for the single differential rate dσ pp /dx F used in [19], where the variables x F and x ′ are defined by In Ref. [19] no information on p ⊥ -distributions was used, and p ⊥ was replaced by some typical value, assumed to bep ⊥ = 1 GeV.
In the present study, the single differential cross section dσ pp /dy can be obtained by integrating (2.11), where the latter approximation arises from µ(p ⊥ ) decreasing much faster than ν(y, p ⊥ ) with p ⊥ . This can be checked for the values of the parameters p 0 , m, n (see Table 1) and within the intervals in p ⊥ and y considered in our study (see Sections 3 and 4). Comparing (A.4) and (A.5) we see that the parametrization (2.11) is consistent with that for the single differential rate (A.4) (or equivalently (A.1)) used in [19]. Let us remark that in (A.5) the typicalp ⊥ may be defined as Using the values of p 0 and m given in Table 1 for J/ψ, we find thatp ⊥ varies in the rangep ⊥ (p 0 , m) = 1.3 − 2.4 GeV, somewhat above the ad hoc valuep ⊥ = 1 GeV used in [19]. However, takingp ⊥ ≃ 2 GeV instead of 1 GeV affects only slightly the value of M ⊥ (M ⊥ ≃ 3.6 GeV instead of 3.2 GeV), with no sizeable effect on the predictions presented in Ref. [19].

B Effective path length vs. centrality class
In this Appendix we explain how the numbers of Table 2 in Section 2.5 have been obtained.

B.1 Number of participants and binary collisions for a given centrality class p-A collisions (LHC case)
In the Glauber model, at a given impact parameter b the number of participating nucleons N p in the target nucleus follows a binomial distribution, Here is the probability for an inelastic collision between the projectile proton and a nucleon of the target nucleus, and σ NN in is the inelastic p-p cross section which we take to be 70 mb at LHC energies. The target nucleus optical thickness T A (b) is normalized as db T A (b) = 1. The denominator in (B.1) is the p-A total inelastic cross section σ pA in and ensures the correct normalization P (N p > 0) = 1 for the total interaction probability.
The number of binary collisions coincides with the number of participating nucleons of the target.

d-A collisions (RHIC case)
The generalization of Eq. (B.1) to the d-A case is straightforward, and is obtained by introducing the distribution P d (r) for the p-n transverse separation r in deuterium, and replacing p pA (b) by the interaction probability p dA (b, r) of a target nucleon with either nucleon of the deuterium projectile, given by Here is the collision probability of a target nucleon with both nucleons of the deuterium projectile. The denominator of (B.2) is the d-A total inelastic cross section σ dA in . The probability (B.4) depends on the N-N inelastic collision probability as a function of the impact parameter p(s). For that profile at c.m. energy √ s = 200 GeV we take the Regge-inspired parametrization p(b) = 1 − exp(−2N e −b 2 /α ) with α = 1.05 fm 2 and N = 1.1, giving the total p-p inelastic cross section σ NN in = ds p(s) = 4.2 fm 2 = 42 mb. The distribution P d (r) is evaluated by assuming a Hulthen form for the deuterium wave function. The thickness functions of all target nuclei considered in our study (Be, Fe, W, Au, Pb) are extracted from low-energy electron-proton scattering experiments [31].
Summing up the probabilities (B.1) and (B.2) we find the threshold values N min p and N max p which saturate approximately 20% of the total interaction probability P (N p ≥ 1) = 1, as described in Section 2.5.
In the d-A case, the number of binary collisions differs from N p by the number of target nucleons which undergo collisions with both nucleons of the deuterium projectile: In the case of triggering on J/ψ production (or any other hard process with a small cross section) the average number of collisions for the participating nucleon gets modified. The hard production process can occur in each of the inelastic collisions of the projectile nucleon with a probability σ J/ψ /σ NN in , 12 so that the hard process probability is P (J/ψ|N p ) = N p σ J/ψ /σ NN in and the joint probability for the hard production and N p participating nucleons in the target is The normalized probability distribution for the number of collisions of the projectile proton in the events tagged by both centrality and J/ψ production is and is independent of the hard process cross section. The corresponding average of N p is The interpretation of (B.8) is straightforward. Unity stands for the target nucleon which participated to the hard process. The second term corresponds to the target nucleons which also undergo an inelastic collision with the projectile, and may thus contribute to the transverse momentum broadening of the cc pair, with the probability σ broad /σ NN in . The effective path length in the target nucleus for the centrality class [N 1 , N 2 ] thus reads where L p is the corresponding length in a proton target. For the minimum bias case doing the summations is trivial and one recovers the expression used in Ref. [19]. In the numerical applications we take L p = 1.5 fm and ρ 0 = 0.17 fm −3 .

d-A (RHIC case)
Compared to the p-A case some complications arise because of the deuterium projectile. What is relevant for the broadening is the number of collisions N t suffered by the nucleon of the deuterium participating to the hard process. This number is not equal to the number of participants in a given event. In the binomial expansion of the first multiplier in the numerator of (B. Correspondingly the number of inelastic rescatterings for the nucleon of the deuterium projectile tagged by the hard process is Interpretation of (B.12) goes along the same line as for (B.8).
The path length of the cc pair for a given centrality class is thus (after substituting p pA (b) = σ NN in T A (b)): 13) For the minimum bias case (N 1 = 1, N 2 = A) the summations in (B.13) are explicit and one again recovers the expression (3.18) of Ref. [19]. The path length for the different centrality classes is given in Table 2.