Nuclear $p_\perp$-broadening of Drell-Yan and quarkonium production from SPS to LHC

The nuclear $p_\perp$-broadening of Drell-Yan and quarkonium ($J/\psi$, $\Upsilon$) production in $\pi$A and pA collisions is investigated. The world data follow a simple scaling from SPS to LHC energies, once the process-dependent color factors are properly taken into account, which allows for the extraction of the transport coefficient in cold nuclear matter. We find that $\hat{q}(x)\propto x^{-\alpha}$ with $\alpha=0.25$-$0.30$. The magnitude of the transport coefficient at $x=10^{-2}$ is $\hat{q}_0=0.051$ GeV$^2$/fm and $\hat{q}_0=0.075$ GeV$^2$/fm, whether $Q\bar{Q}$ pairs are assumed to be produced as color octet or color singlet states, respectively. The relation between nuclear broadening data and the (CT14) gluon density is also investigated.


Introduction
The multiple scattering incurred by energetic quarks and gluons propagating in a QCD medium is responsible for both induced gluon radiation -leading to parton energy loss -and transverse momentum broadening. While the former is responsible for the jet quenching phenomena discovered in heavy ion collisions at RHIC and at LHC (see Refs. [1][2][3][4] for reviews), observing the latter is more delicate. Measuring dijet azimuthal correlations in proton-nucleus [5,6] and nucleus-nucleus collisions could in principle allow for probing transverse momentum broadening, respectively in cold nuclear matter and hot quark-gluon plasma [7][8][9][10]. In practice, however, these correlations may prove more sensitive to medium-independent Sudakov radiation rather than to multiple scattering in a medium and thus deserve detailed scrutiny [11].
A perhaps simpler observable, namely the nuclear broadening of single inclusive particle production in hA collisions, offers a direct access to the saturation scale, hence a measure of multiple scattering in nuclei [12][13][14]. In this article, we perform a systematic study of Drell-Yan (DY) and quarkonium (J/ψ, Υ) nuclear broadening data in πA and pA collisions, from fixed-target (SPS, FNAL) to collider (RHIC, LHC) energies. It is found that all measurements, once scaled with appropriate color factors, exhibit a simple dependence expected from perturbative QCD. This data-driven analysis allows in turn for the determination of the transport coefficient in cold nuclear matter, or equivalently the saturation scale in large nuclei and its energy dependence. The model for the nuclear broadening is presented in section 2 before results are shown and discussed in section 3. Conclusions are drawn in section 4.
2 Model and data 2.1 Transverse momentum broadening

Asymptotic parton
The accumulated transverse momentum acquired by an asymptotic parton (in color representation R) traversing a medium of length L is given by [12] ∆p 2 ⊥ =q R L , (2.1) whereq R is the transport coefficient of that medium. 1 It is given by µ 2 /λ R where µ is the typical transverse momentum exchange between the parton and each scattering center, and λ R is the parton mean free path in the medium. Definingq as the gluon transport coefficient and using C R λ R = N c λ g , the broadening becomes where C R is the color charge of parton R (C F = (N 2 c − 1)/2N c for a quark, C A = N c for a gluon).

Hard QCD process on a nuclear target
Let us now consider the more realistic case of the production of a single massive particle produced by a hard QCD process in a hadron-nucleus collision. 2 When the coherence length of the hard process is small, typically smaller than the internucleonic distance, coh 1 fm, the process as seen in the nuclear rest frame can be sketched as in Figure 1: (i) an incoming parton (in color representation R) stemming from the incoming hadron experiences transverse momentum broadening over a distance z, at which it participates to the hard process; (ii) an outgoing particle (in color representation R ) is produced on a short length scale coh and propagates in the medium over a length L − z, until it exits the nucleus. In this picture, the nuclear transverse momentum broadening in hA collisions compared to hp collisions is simply given by 3 where the color factor C is given by the half sum of the Casimir factors of the initial and final-state particle, 4 Assuming a hard sphere nuclear density of radius R A = 1.12 × A 1/3 fm, the average medium length in a large nucleus is given by L A = 3/2 R A after averaging uniformly over the position of the production process. The length scale in a proton target is set to L p = 1.5 fm, as in Refs. [18,19].
In the limit of large coherence length, coh L, the broadening given by Eq. (2.3) with the generic color rule Eq. (2.4) may not hold for all hard QCD processes. In that limit, the outgoing particle (in the present context, a virtual photon or a compact QQ pair) detected experimentally would occur, at the amplitude level, from the radiation of the incoming parton either long before or long after crossing the nuclear target in its rest frame. In the case of initial-state (respectively, final-state) radiation in both the amplitude and complex conjugate amplitude, the broadening would depend solely on the Casimir of the outgoing particle C R (respectively, that of the initial parton C R ), while the broadening associated to the interfering terms between initial and final state radiation would have no simple expression in terms of Casimir factors. Consequently, Eq. (2.3) may not be appropriate to describe the broadening of Drell-Yan lepton pairs at high energy, as the initial parton (a quark) and the outgoing particle (the virtual photon) carry different color charges. In any case, measurements of DY broadening are not available yet at large coherence length, i.e. at RHIC and LHC, despite the preliminary results being reported by the PHENIX experiment [20] which are not included in this analysis. On the contrary, when the incoming and outgoing particles carry the same color charge, as would be the case for an octet QQ pair induced by an incoming gluon (C R = C R ), Eq. (2.3) is likely to hold (i.e., C = N c ). We shall therefore assume in the following that the broadening of quarkonia at high energy (RHIC, LHC) is given by Eq. (2.3).

Process dependence
In this section we discuss the specific color factor C expected for Drell-Yan and quarkonium production in πA and pA collisions.

Drell-Yan
At Born level, DY lepton pairs are produced in hadron-nucleus collisions by the annihilation of a quark (respectively, an antiquark) from the projectile hadron with an antiquark (respectively, a quark) from the nuclear target. Since the virtual photon (or the lepton pair) does not experience multiple scattering in the target, C R = 0, the expected color factor is thus C DY = C F /2, independently of the nature of the projectile hadron (here, a pion or a proton). At order O (α s ), the QCD Compton process, qg → qγ , may come into play. However, it is not expected to dominate the cross section at low p ⊥ M DY , where M DY is the invariant dilepton mass. Moreover, at forward dilepton rapidity (as is the case for the DY data presently analyzed), this process is dominated by the fusion of an incoming quark with a gluon from the target (as q(x 1 )g(x 2 ) g(x 1 )q(x 2 )) leading to the same color factor as for the qq annihilation process.

Quarkonium
Which color factor to use for quarkonium production is more delicate for two reasons. At leading order in α s , both gluon fusion (gg → QQ) and quark-antiquark annihilation (qq → QQ) processes are expected to contribute to heavy-quark pair production, which would lead, respectively, to C R = N c and C R = C F in Eq. (2.4). In pA collisions, however, quarkonium production should be dominated by gluon fusion unless the longitudinal momentum fraction becomes very large, x F 0.5. This statement depends mostly on the momentum distributions of quarks and gluons, and appears to be true in non-relativistic QCD (NRQCD) or in the Color Evaporation Model (CEM) [21], leading to C R = N c in (2.4). The qq annihilation process is more significant in πA collisions as the pion carries valence antiquarks. Using the pion Table 1. Color factors assumed in the present analysis for Drell-Yan and color octet quarkonium production in πA and pA collisions. The assumption regarding the case of color singlet quarkonium production is discussed in section 3.3.
PDF recently extracted in Ref. [22] in the CEM at leading order, we find that the qq annihilation channel actually dominates the inclusive J/ψ production in πA collisions (hence, C R = C F ), at least in the x F range considered in this analysis. In addition, the quarkonium production process is still poorly understood. In particular it is not clear on which length scale ( octet ) the color octet heavy-quark pair neutralizes its color, either octet ∼ coh (as in the Color Singlet Model) or octet coh (NRQCD, CEM). In this study we shall assume that the QQ pair remains in a color octet state during its entire propagation throughout the nucleus, leading to C R = N c . For completeness, the assumption of color singlet QQ production is also discussed in section 3.3. The color factors assumed in this study for DY and quarkonium production in hadron-nucleus collisions are summarized in Table 1.

Transport coefficient
Apart from the Casimir scaling properties discussed in the previous section, the other crucial ingredient which governs the nuclear broadening of hard QCD processes iŝ q, the gluon transport coefficient in cold nuclear matter. It is related to the gluon distribution xG(x, Q 2 ) inside each nucleon of the target [12], where ρ is the nuclear density.
In the case of a hard parton produced inside the medium (consider for instance a quark produced in deep inelastic scattering at large Bjorken-x), typically when the coherence length is small, coh L A , the value of x at which the gluon distribution needs to be probed is x ∼ x A ≡ 1/(2m N L A ) > 10 −2 [12]. On the contrary, when the hard parton is produced by a QCD process coherent over the whole nucleus, coh L A , x should be given by x ∼ 1/(2m N coh ) = x 2 , where x 2 is the momentum fraction carried by the parton probed in the nuclear target [19]. The scale Q 2 which controls both the running coupling α s and the QCD evolution of xG should be of the order of the broadening itself, Q 2 ∼ ∆p 2 ⊥ , yet in practice this semi-hard scale should be frozen when ∆p 2 ⊥ becomes too small, ∆p 2 ⊥ 1-2 GeV 2 . At small x 10 −2 , the gluon distribution exhibits a power law behavior, xG(x, Q 2 ) ∝ x −α , where the exponent α 0.2-0.3 depends only slightly on the resolution scale Q. Neglecting the Q 2 dependence ofq in Eq. (2.5), the transport coefficient can thus be modelled as [19] whereq 0 is the transport coefficient at x = 10 −2 . Note that at large x 2 > x A ,q acquires a mild dependence on the size of the nucleus from the L dependence of x A . Casting (2.6) in (2.3) leads to At small x 2 < x A , L p coincides with L p , while at large x 2 > x p the factor accounts for the fact that the transport coefficient should be evaluated at two different values of x in the proton and in the nucleus targets. The value of x 2 is given by In the present approach, the absolute magnitudeq 0 and the slope α are the only parameters which need to be extracted from the measurements. The relationship (2.5) between the nuclear broadening and more realistic gluon distribution functions will be investigated in section 3.2.

Other nuclear effects
It has been assumed so far that parton multiple scattering in nuclei is the only process which leads to finite nuclear p ⊥ -broadening, ∆p 2 ⊥ = 0. However, other cold nuclear matter effects such as fully coherent energy loss (FCEL) [23][24][25] or nuclear parton distribution functions (nPDF) [26][27][28][29][30] may distort the shape of particle p ⊥ -spectra in hadron-nucleus collisions with respect to hadron-proton collisions.
The p ⊥ dependence of the nuclear production ratio, due to fully coherent energy loss has been investigated in the case of quarkonium [31] and light hadron [32,33] production. The average fully coherent energy loss is suppressed by one power of the hard scale, ∆E ∝ 1/p ⊥ at p ⊥ M , making R pA (p ⊥ ) a growing function of p ⊥ in the rapidity intervals presently discussed. Consequently, the quarkonium p ⊥ -spectra prove harder in hA collisions than in hp collisions, thus leading to a positive 'nuclear broadening', ∆p 2 ⊥ FCEL > 0. Unlike quarkonium production, no FCEL is expected in the Drell-Yan channel [34], which may nonetheless be sensitive to initial-state energy loss. These effects, however, prove tiny either when x F is not too large or at high collision energies [35].
Let us now discuss the role of nPDF effects on nuclear p ⊥ -broadening. At large p ⊥ M , the momentum fraction probed in the nuclear target depends on the transverse momentum of the detected particle, x 2 = x 2 (p ⊥ ). Assuming for simplicity that only one parton species from the nucleus (call it flavour i) dominates the cross section, the p ⊥ spectrum in hA collisions may be given by (neglecting here broadening effects) is the nPDF ratio of parton flavour i in the nucleus A over that in a proton, and the factorization scale Q might itself depends on p ⊥ . When R A i varies little with p ⊥ compared to the weighted spectrum p 2 ⊥ dσ hp /dp ⊥ , the shape of the spectrum is not affected by nPDF corrections but only its magnitude. As a consequence, no net effect of nuclear parton densities on the broadening is expected. This may be the case in the deep shadowing region at very small x 2 , typically x 2 10 −4 , and in the vicinity of the anti-shadowing region, 0.05 x 2 0.2. 6 In between these two domains, R i is often a growing function of x 2 (and thus of p ⊥ ), making the particle spectrum harder in hA collisions compared to hp collisions. As a consequence, the sole nPDF effects would lead to a positive contribution to the nuclear broadening, ∆p 2 ⊥ nPDF > 0, if the typical transverse momenta contributing to p 2 ⊥ correspond to these domains in x 2 . Conversely, at larger x 2 0.2, the possible decrease of R A i with x 2 due to the EMC effect would lead to softer spectra in hA collisions, leading to a negative contribution to the nuclear broadening, ∆p 2 ⊥ nPDF < 0. At fixed-target collision energies, for which p ⊥ M , the momentum fraction depends only mildly on p ⊥ , therefore no strong nPDF effects on the broadening are expected.
Quantitative results on the FCEL and nPDF effects to the nuclear broadening are discussed in Appendix A. These prove significantly less than in data, which confirms that the measured nuclear broadening can be mostly attributed to multiple scattering effects, rather than to FCEL or nPDF.

Data
In the present analysis, all available data on the nuclear broadening of Drell-Yan, J/ψ, and Υ production have been used  The analyzed data sets are summarized in Table 2. At SPS (NA3, NA10), measurements have been performed in πA and pA collisions, which allows for investigating the color charge dependence of J/ψ nuclear broadening; see Table 1. The E772 results performed in pA collisions on different nuclear targets are sensitive to the path-length dependence of both DY and Υ nuclear broadening. Finally, the RHIC and LHC measurements are carried out on a single nuclear target (Au and Pb, respectively) but in different rapidity ranges, thus probing the x 2 dependence of the transport coefficient, as emphasized e.g. in Ref. [14]. In particular, the forward measurements at RHIC (1.2 < y < 2.2) and at LHC (2 < y < 4, for the LHCb experiment) correspond to data taken at the smallest x 2 values, x 2 ≈ 3 × 10 −3 and x 2 ≈ 2 × 10 −5 , respectively.
At the LHC, the J/ψ nuclear p ⊥ -broadening in minimum bias pA collisions have not been released by the ALICE and LHCb collaborations. The extraction of ∆p 2 ⊥ from the measurements of the absolute p ⊥ -spectra, as well as a new extraction from PHENIX data [39], are detailed in Appendix B.

Scaling
In order to test the above picture, the world data on nuclear transverse momentum broadening of Drell-Yan, J/ψ and Υ production is plotted in Figure 2 as a function of C ∆L/x α , see Eq. (2.7). Remarkably, the measurements performed in pA and πA collisions, for different particle species and on a large range of collision energies from   Figure 2. Scaling of the nuclear p ⊥ -broadening in DY and quarkonium production using the color assumptions in Table 1.
SPS to LHC, exhibit clearly the scaling property. As expected, the lowest values of ∆p 2 ⊥ are obtained in the DY process (C = C F /2) at low collision energy (hence, larger x 2 for whichq is lower). On the contrary, the largest broadening is observed for quarkonium production at the LHC, because of both the large color factor (C = N c ) and the small-x 2 rise of the gluon distribution,q(x 2 ) ∝ x −α 2 at small x 2 .
Let us now discuss the values of the two parameters of the model obtained from the fit to the data. The best fit (χ 2 /ndf = 1.3) leads to the small-x 2 exponent of the transport coefficient α = 0.25 ± 0.01, as can be seen from the χ 2 /ndf profile shown as a dotted line in Fig. 3. This value of α proves in very good agreement with the 'geometrical scaling' fits of HERA data [42], which is quite remarkable given the different observables. It is also nicely compatible with the increase of the saturation scale deduced from elliptic flow measurements in heavy ion collisions from RHIC to LHC [43]. Recently a study based on the higher-twist framework, aiming at the extraction of the transport coefficient from a global fit of semi-inclusive DIS and pA collisions data led to a slightly smaller exponent,q(x 2 ) ∝ x −0.17 2 [44].
The absolute value of the transport coefficient obtained from the fit,q 0 = 0.051± 0.002 GeV 2 /fm, coincides with the perturbative estimate by BDMPS in Ref. [12]. It  is also consistent with estimates from fully coherent energy loss effects [19]. 7 The quoted uncertainties on the two parameters follow from the χ 2 analysis with the strict tolerance criterion of ∆χ 2 = +1. However, it is clear that these uncertainties should be taken as lower estimates, as other sources of 'systematic' uncertainties inherent in the model are likely to affect these results. One such assumption is the value of x at which the transport coefficient is frozen, x A = 1/(2m N L A ), see (2.6). Yet motivated on physical grounds by the uncertainty principle, this expression gives at most the magnitude for x A . Varying x A by a factor of two with respect to that estimate leads respectively toq 0 = 0.055±0.002 GeV 2 /fm andq 0 = 0.044±0.001 GeV 2 /fm, still with a rather good agreement with data (χ 2 /ndf = 1.9 and χ 2 /ndf = 1.4, respectively), thus leading to an additional uncertainty onq 0 of about 10%. It has also been checked that evaluatingq(x 2 ) at different values of x 2 , e.g. x 2 = M 2 + p 2 ⊥ 1/2 / √ s × e −y at various p 2 ⊥ , affects only marginally (by less than 5%) our results. The largest source of uncertainty onq 0 is related to the assumption on the color state of the propagating QQ state; it is discussed more specifically in section 3.3. It has been pointed out that the broadening may have a different expression at large coherence length for specific QCD processes, typically when the color charge of the incoming parton and that of the outgoing particle differ (e.g. Drell-Yan or color singlet QQ production) while (2.3) should hold at any coherence length in the case of color octet QQ production, as assumed here. Nevertheless, in order to check the consistency of our approach at low and at high energy, a fit to a subset of data at small coherence length has been performed, taking as a criterion coh < L/2 which excludes the LHC J/ψ data and the RHIC J/ψ measurements at mid and forward rapidity. The fit leads toq 0 = 0.050±0.002 GeV 2 /fm (χ 2 /ndf = 0.8 with α = 0.25), a value which is nicely consistent with results from the full data sets, showing that the RHIC and LHC data do not alter the value extracted from low energy measurements.
High energy data are nevertheless key in order to extract precisely α, the small-x 2 exponent of the transport coefficient. Indeed, at low energy the transport coefficient (2.6) is independent of x 2 and solely depends on x A ≡ 1/(2m N L). 8 Since most of the nuclei involved in this analysis have similar atomic masses hence similar medium lengths (A 200 corresponding to x A 10 −2 ), the transport coefficient used in the low-energy fit is a constant independent of α,q(x) q 0 . These low-energy data alone are thus unable to constrain the parameter α, as can be seen in Figure 3 (solid line) which shows a rather flat χ 2 /ndf profile as a function α. The results of the fits, either to the full data set or to the smallcoh subset of data, are summarized in Table 3.  Table 3. Results of the fits to all and selected ( c < L/2) data sets, assuming color octet or color singlet QQ production.
On top of the global fits, the values ofq 0 have been extracted from each data point using (2.3) and (2.6). Within each experiment, the different values correspond either to different nuclei (E772), different collision energies (NA3, NA10), and to different rapidity bins (PHENIX, ALICE, LHCb), see Table 2. Results plotted in Figure 4 show a remarkable consistency (as could be anticipated from Figure 2), pointing to a common value for the transport coefficientq 0 . The weighted average ofq 0 is also consistent with the estimates above. It is interesting to note that the   Drell-Yan measurements by NA10 lie slightly above the average. The largest tension is observed with the LHCb measurement at backward rapidity (−4 < y < −2), which might be partly attributed to FCEL or nPDF effects, as discussed in Appendix A.

Relationship with the gluon distribution
As discussed in section 2.3, the transport coefficient is directly proportional to the gluon distribution of the nucleus, Eq. (2.5), which led to Eq. (2.6), tested in the previous section. The large range of variation in the value of x, from x 2 × 10 −2 at low collision energy to x 2 × 10 −5 at the LHC at forward rapidity, 9 makes it possible to check further the relationship betweenq and xG(x) by comparing data to actual proton PDF sets. From (2.5), the broadening reads where the proportionality constant is related toq 0 which is now the only parameter, and x ≡ min(x p , x 2 ). As mentioned in section 3.2, the scale Q 2 appearing in (3.1) is expected to be of the order of the broadening itself, Q 2 = λ ∆p 2 ⊥ with λ = O (1). At such semi-hard scales (look at Figure 2, in which most data points lie below ∆p 2 ⊥ 1 GeV 2 ), parton densities are not defined and Q 2 should be frozen at the input scale, Q 2 = min(Q 2 0 , λ∆p 2 ⊥ ) defined in the global fit analysis (typically, Q 2 0 ∼ 1-2 GeV 2 ). The fit has been performed using the leading-order 10 CT14 PDF set [45]. A bit surprisingly, the choice λ = 1 leads to a rather poor agreement with data. Due to the quasi-absence of QCD evolution from Q 2 0 to ∆p 2 ⊥ , the gluon distribution does not rise sufficiently at small x to account for the measurements. Better agreement is nevertheless found for larger values of λ, e.g. λ 2.5 (χ 2 /ndf = 1.2) as can be seen in Figure 5 which compares data and the assumption (3.1). Moreover the obtained value for the transport coefficient at x = 10 −2 ,q 0 = 0.049 ± 0.002 GeV 2 /fm, is compatible with our previous estimate based on a simpler model for the gluon distribution.

Color singlet production
The present lack of understanding regarding the production process of quarkonium states makes it difficult to predict the nuclear broadening of J/ψ and Υ states. It has been assumed in this analysis that the QQ pair neutralizes its color on a long timescale, as in the CEM or NRQCD, and for which Eq. (2.3) appears appropriate at both low and high energies. In this section, we shall assume that QQ states turn color singlet on a short timescale, in the spirit of the color singlet model (CSM). In the CSM, the leading order process for quarkonium production is gluon fusion, thus leading to the color factor C = N c /2 in both pA and πA collisions 11 (see Table 1).
As Eq. (2.3) may not hold at large coherence length in the case of color singlet QQ production, we first compare the model expectations to the DY and quarkonium data at small coherence length, coh < L/2. The fit of this subset of data leads to a transport coefficientq 0 = 0.074 ± 0.003 GeV 2 /fm (χ 2 /ndf = 1.3 at α = 0.30). Without much surprise, the transport coefficient extracted assuming color singlet QQ production proves somewhat larger than assuming color octet QQ pairs, in order to compensate for the smaller color factor, see Table 1. As in the case of color octet QQ production, these low energy data alone do not allow for the determination of the small-x 2 exponent α; see the flattish χ 2 /ndf profile in Fig. 3 (dashed line). For completeness, a fit of the full data set has been performed under the color singlet assumption. The fits leads toq 0 = 0.075 ± 0.003 GeV 2 /fm (however with a rather large χ 2 /ndf = 2.0), which is also consistent with the value obtained from low energy data. The use of the whole data set now allow for a precise extraction of α, α = 0.30 ± 0.02 (Figure 3, dash-dotted line). Interestingly, this estimate is comparable in magnitude with our earlier result assuming color octet QQ states, and also compatible with small-x deep inelastic scattering data [42]. This result should however be taken with care as Eq. (2.3) may not be valid at high energy in the case of color singlet production.

Conclusion
The transverse momentum nuclear broadening of Drell-Yan lepton pair and quarkonium production has been investigated systematically, for different nuclear targets and different systems (pA and πA collisions), from SPS to LHC energy. Within a model that includes respectively the dependence on the medium length, the processdependent color factors, and the x dependence of the transport coefficient (or equivalently that of the saturation scale), a simple scaling is expected and is observed in the data.
This allows for the determination of the transport coefficient of cold nuclear matter at x = 10 −2 ,q 0 = 0.051 GeV 2 /fm, with a systematic uncertainty estimated to be 10%. Moreover, the best fit to data points to the x-dependence,q(x) ∼ x −0.25 . This would correspond to a saturation scale in a large nucleus (either Au or Pb) of Q s = 0.7 GeV at RHIC at mid-rapidity and Q s = 1.1 GeV at LHC at mid-rapidity (using here x 2 M J/ψ / √ s). The relation expected in perturbative QCD between the transport coefficient and the gluon distribution in a proton is explored further using CT14 LO gluon distribution. Good agreement is reported provided the hard scale entering xG(x, Q 2 ) and α s (Q 2 ) is Q 2 2.5 ∆p 2 ⊥ , that is slightly above the BDMPS prescription Q 2 = ∆p 2 ⊥ [12]. The analysis has also been carried out assuming that QQ pairs turn color singlet on a short timescale. In this case, the full dataset would favor a naturally larger transport coefficient,q 0 = 0.075 GeV 2 /fm, albeit with a similar small-x exponent, α = 0.3.
The present picture could be further tested with additional measurements. At the LHC, Drell-Yan data at backward/forward rapidity by LHCb and at midrapidity by ATLAS and CMS would be extremely valuable. These data would shed light on the DY broadening at small values of x 2 , that is when the coherence length for the hard process is significantly larger than the medium length (the RHIC forward DY data may also be important in this respect [20]). In addition, at lower energy, the simultaneous measurement of DY and J/ψ nuclear broadening in πA collisions at forward rapidity by the COMPASS experiment [46], or in pA collisions at backward rapidity by the LHCb-SMOG experiment [47,48], should also bring essential information on the value of the transport coefficient as well as on quarkonium formation dynamics.
A FCEL and nPDF effects on ∆p 2 ⊥ The effects of FCEL and nPDF on the shape of J/ψ p ⊥ -spectra in pA collisions, hence on the nuclear broadening ∆p 2 ⊥ , are investigated. Neglecting multiple scattering and focusing on the sole FCEL/nPDF effects, the cross section in pA collisions can be modelled as (see section 2.4), The FCEL quarkonium nuclear production ratio R FCEL pA is computed from Ref. [19], while the gluon nPDF ratio R A g is given by EPPS16 [29]. The choice made here for x 2 and Q that enter R A g interpolates between the 2 → 1 and 2 → 2 kinematics expected at low and high p ⊥ . The double differential J/ψ production cross section in pp collisions is parametrized as where the values of the parameters N , n, m and p 0 were extracted from √ s = 200 GeV and √ s = 7 TeV data in Ref. [31]. Using the spectrum (A.2) in (A.1) allows for computing ∆p 2 ⊥ FCEL and ∆p 2 ⊥ nPDF . Calculations are performed at RHIC and at LHC ( √ s = 8.16 TeV), in the rapidity acceptance of the PHENIX experiment (1.2 < |y| < 2.2 and |y| < 0.35) and the LHCb experiment (−4.5 < y < −2.5 and 2 < y < 4), respectively. 12 Results are summarized in Table 4. At RHIC, the values of ∆p 2 ⊥ due to FCEL are small, ∆p 2 ⊥ FCEL = 0.1 GeV 2 in all rapidity bins. At the LHC, FCEL alone leads to transverse momentum broadening ranging from ∆p 2 ⊥ = 0.2 GeV 2 at backward rapidity (y = −3.5) to ∆p 2 ⊥ = 0.5-0.6 GeV 2 at forward rapidity (y = 3), consistent with the fact that FCEL effects are more pronounced in the proton fragmentation region. The quoted uncertainties arise from the variation of the transport coefficient fromq 0 = 0.050 toq 0 = 0.075 GeV 2 /fm, consistently with the results obtained in this analysis.
Let us now comment on the nPDF effects. At RHIC, the nPDF contribution to ∆p 2 ⊥ at backward rapidity is slightly negative as the region between anti-shadowing and the EMC effect softens p ⊥ -spectra in pA collisions with respect to pp collisions (see discussion in section 2.4). The quoted uncertainty is determined after the variation of all EPPS16 nuclear member sets. At mid-rapidity and at forward rapidity, ∆p 2 ⊥ nPDF is positive in the two rapidity bins (0.1 < ∆p 2 ⊥ nPDF < 0.4 GeV 2 ). At LHC, the values of ∆p 2 ⊥ range from ∆p 2 ⊥ = 0.1 to ∆p 2 ⊥ = 0.6 GeV 2 at backward rapidity and from ∆p 2 ⊥ = 0.1 to ∆p 2 ⊥ = 0.7 GeV 2 at forward rapidity, with a significant correlation observed for each EPPS16 member set between the two rapidity intervals (correlation coefficient of 0.7).

Experiment
System y range ∆p 2 ⊥ FCEL (GeV 2 ) ∆p 2 ⊥ nPDF (GeV 2 ) PHENIX dAu −2.2 < y < −1. Although they might a play a role on the values of ∆p 2 ⊥ , this study reveals that both FCEL and nPDF effects appear to be small (yet with a large uncertainty in the latter case) compared to the values observed in data, suggesting that multiple scattering is the leading effect, as assumed throughout the present analysis.

B Extraction of ∆p 2 ⊥ at RHIC and LHC
In practice, the average transverse momentum squared defined as cannot be extracted from data. Instead, p 2 ⊥ is often determined from the sum over the experimental bins (with lower and upper edges, p i−1 ⊥ and p i ⊥ ), taken as an approximation of Eq. (B.1). This latter expression has however two drawbacks. When the bins are too wide, as is often the case at large p ⊥ , the typical p i ⊥ value at which the transverse momentum should be evaluated in the bin [p i−1 ⊥ , p i ⊥ ], for instance the average or the median in that bin, may significantly affect the value of p 2 ⊥ bins . In addition, Eq. (B.2) should be used only when the highest p ⊥ value reached in the experiment, p max ⊥ = p N bins ⊥ , is large enough so that p 2 ⊥ bins is independent of this upper cutoff, within the experimental uncertainty.
Whenever available, we use in the present article the values of ∆p 2 ⊥ published in the experimental analyses. This is the case of all fixed-target experiments. The specific case of PHENIX data is specifically discussed hereafter. At the LHC the values of ∆p 2 ⊥ in minimum bias pPb collisions have not been published, neither by ALICE nor by LHCb. The extraction of ∆p 2 ⊥ from these data is discussed below.
PHENIX -The PHENIX experiment published p 2 ⊥ values for J/ψ production in pp and dAu collisions at √ s = 200 GeV in different rapidity bins [39]. Using the absolute cross sections measured in pp and dAu collisions [39,49], we have checked that the broadening values quoted in [39] can be recovered when using the median transverse momentum,p i ⊥ = (p i−1 ⊥ +p i ⊥ )/2, in Eq. (B.2). However, because of the wide experimental bins at large p ⊥ , the values of p 2 ⊥ bins in pp and dAu collisions, hence the nuclear broadening ∆p 2 ⊥ , is significantly affected when replacing the median by the average transverse momentum in each bin 13 in Eq. (B.2). In order to circumvent the drawbacks of using Eq. (B.2), data have been fitted using the so-called Kaplan distribution, dσ fit (p ⊥ ) p ⊥ dp ⊥ = N p 2 0 p 2 0 + p 2 ⊥ m (B.3) 13 This value is not given in [39] but was estimated using the fit (B.3) in each bin. named after Ref. [50]. Once the parameters m and p 0 are known (the parameter N is irrelevant when computing p 2 ⊥ ), the value of ∆p 2 ⊥ can be determined using the proper definition Eq. (B.1).
LHCb -The same procedure, namely using (B.3) in (B.2), is applied to J/ψ production cross sections measured by LHCb in pp and pPb collisions at √ s = 8.16 TeV [41].

ALICE -
The ALICE experiment has released the values of ∆p 2 ⊥ for J/ψ production in pPb (and pp) collisions at √ s = 5.02 TeV [40] in different centrality bins. Here the weighted average over the N centrality bins has been carried out, where σ C are the measured p ⊥ -integrated J/ψ production cross section in each centrality bin. The values of ∆p 2 ⊥ extracted from PHENIX, LHCb and ALICE are given in Table 5.