Exclusive $J/\psi$ photoproduction in ultraperipheral Pb+Pb collisions at the LHC to next-to-leading order perturbative QCD

We present the first next-to-leading-order (NLO) perturbative QCD (pQCD) study of rapidity-differential cross sections of coherent exclusive photoproduction of $J/\psi$ mesons in heavy-ion ultraperipheral collisions (UPCs) at the LHC, $d\sigma/dy(\text{Pb} + \text{Pb} \rightarrow \text{Pb} + J/\psi + \text{Pb})$. For this, we account for the photon-nucleon NLO cross sections at the forward limit, the $t$ dependence using a standard nuclear form factor, and the photon fluxes of the colliding nuclei. Approximating the generalized parton distributions with their forward-limit parton distribution functions (PDFs), we quantify the NLO contributions in the cross sections, show that the real part of the amplitude and quark-PDF contributions must not be neglected, quantify the uncertainties arising from the scale-choice and PDFs, and compare our results with ALICE, CMS and LHCb $J/\psi$ photoproduction data in Pb+Pb UPCs, exclusive $J/\psi$ photoproduction data from HERA, and LHCb data in p+p. The scale dependence in $d\sigma/dy(\text{Pb} + \text{Pb} \rightarrow \text{Pb} + J/\psi + \text{Pb})$ is significant, but we can find a scale-choice that reproduces the Pb+Pb UPC data both at 2.76 and 5.02 TeV collision energies. This process has traditionally been suggested to be a direct probe of nuclear gluon distributions. We show that the situation changes rather dramatically from LO to NLO: the NLO cross sections reflect the nuclear effects of both gluons and quarks in a complicated manner where the relative signs of the LO and NLO terms in the amplitude play a significant role.


I. INTRODUCTION
Ultraperipheral collisions (UPCs) are collisions of hadrons or nuclei which take place at large impact parameters in such a way that only the electromagnetic field of one of the colliding particles interacts with the other particle [1][2][3]. Coherent photoproduction of J/ψ heavy vector-mesons in UPCs of Lead nuclei at the CERN Large Hadron Collider (LHC), the exclusive process Pb+Pb → Pb+J/ψ+Pb, has been suggested to be an efficient direct probe of collinear nuclear gluon distributions, g Pb (x, Q 2 ), at factorization scales of the order of the vector-meson mass, Q 2 = O(M 2 V ), and small longitudinal-momentum fractions x = O(M 2 V /W 2 ), where W is the photonnucleon center-of-momentum-system (c.m.s.) energy [4][5][6][7][8][9][10][11]. This exciting possibility derives from the fact that in such an exclusive process of no hadronic activity, one of the colliding nuclei serves as a source of equivalent real Weizsäcker-Williams photons which probe a colorsinglet gluon-or quark-initiated ladder from the other nucleus via formation of a heavy quark-antiquark pair. As first discussed by Ryskin in Ref. [12] in the context of the free-proton process γ + p → J/ψ + p, in the leading order (LO) perturbative QCD (pQCD) only the gluon-ladder processes contribute, and neglecting the longitudinal-momentum imbalance (skewedness) in the ladder and the subleading real part of the amplitude, the forward scattering amplitude factorizes into a calculable hard part and g p (x, Q 2 ). Thus the cross section of J/ψ * topi.m.o.loytainen@jyu.fi becomes proportional to [g p (x, Q 2 )] 2 , making the process a very promising one for probing the gluon distribution. This idea has then been transferred to ultraperipheral nucleus-nucleus collisions (UPCs) in e.g. Refs. [4,5]. Also Monte Carlo event simulations of this process in the UPCs have been developed, such as STARlight [13] and SuperChic [14]. Exclusive photoproduction of J/ψ has also been widely studied in the dipole picture, especially in the high-energy Color-Glass-Condensate approximation of QCD, see e.g. Refs. [15][16][17][18][19][20][21][22][23][24][25].
With the experimental data being released from the LHC, the situation is becoming ever more interesting. Firstly, the exclusive coherent J/ψ photoproduction cross sections involving real photons have been measured in electron-proton collisions at the DESY-HERA collider by the H1 [26] and ZEUS [27] collaborations, and extracted also from the LHCb measurements of the process p + p → p + J/ψ + p at the LHC [28,29]. For detailed NLO pQCD studies of these, see e.g. Refs. [30][31][32][33][34][35]. From the viewpoint of the UPCs, these data sets offer also an importantly long lever arm in the photonproton c.m.s. energy W for cross-checking the pQCD calculations and understanding the necessary modeling input. Secondly, in Pb+Pb UPCs at the LHC, the AL-ICE collaboration has measured the rapidity-differential cross section of Pb + Pb → Pb + J/ψ + Pb both at midrapidity [36,37] and at forward/backward rapidities [38,39] at nucleon-nucleon c.m.s. energies √ s NN = 5.02 TeV and 2.76 TeV. The CMS collaboration has performed the corresponding measurement at √ s NN = 2.76 TeV in one off-central rapidity bin that lies conveniently just between the ALICE rapidity bins [40]. The LHCb collaboration has recently released their 5.02 TeV data at arXiv:2203.11613v1 [hep-ph] 22 Mar 2022 forward/backward rapidities [41], overlapping with the ALICE rapidity region. Very interestingly, however, the ALICE and LHCb forward/backward-rapidity 5.02 TeV data sets do not seem to be fully compatible with each other, which clearly calls for further analyses. Until now, exclusive J/ψ photoproduction in ultraperipheral nuclear collisions has been studied only to LO pQCD. Now that the LHC experiments are measuring these cross sections to an increasing accuracy, and hopefully also for other UPC systems than Pb+Pb in the future [42], it is clearly of high priority to extend the theory calculations to NLO pQCD. In particular we wish to study whether/how this process could be included in the global analyses of nuclear PDFs, such as in Refs. [43][44][45][46][47], in the future. These are the main motivations for our present NLO study. Also interestingly, so far the LO pQCD, or dipole picture, calculations have not been able to reproduce simultaneously the mid-and forward/backward-rapidity data, see e.g. [36]. This, together with the mentioned incompatibility between the LHCb and ALICE data, serves also as further motivation for our current NLO pQCD study.
The NLO pQCD calculation of cross sections for exclusive photoproduction of heavy vector mesons V off the free proton, σ(γ + p → V + p), using collinear factorization at the amplitude level, has been performed first by Ivanov et al. in Ref. [30], followed then by other groups in Refs. [31,35,[48][49][50]. To be exact, collinear factorization here refers to the factorization of the amplitude to calculable NLO pQCD pieces and to the generalized parton distributions (GPDs) [51] which at the forward limit relax into the usual PDFs [52]. If such a limit is not assumed, then the GPDs have to be modeled in some way, e.g. as suggested in Refs. [53][54][55][56][57][58][59]. As is shown already in Ref. [30], the full NLO calculation of coherent exclusive photoproduction of J/ψ mesons in γ + p collisions, which includes both the imaginary and real parts of the amplitude precisely as they are, and assumes a certain model for the gluon-and quark-GPDs [55], depends rather heavily on the choice of the renormalization/factorization scale, Q = O(M J/ψ ), while for the photoproduction of Υ mesons, which probes a higher scale Q = O(M Υ ), the situation improves somewhat. Discussion of a systematic procedure for diminishing the scale dependence in the NLO calculation of exclusive J/ψ photoproduction in γ + p collisions can be found in [32][33][34][35], but in the present exploratory NLO study for the nuclear UPCs we do not follow this avenue.
In the current paper, we present the first NLO pQCD study of exclusive photoproduction of J/ψ mesons in ultraperipheral Pb+Pb collisions at the LHC, with collinear factorization at the amplitude level. Exploiting the analytic results of the impressive calculation of Ref. [30], we have built a numerical code of our own for the rapidity-differential J/ψ photoproduction UPC cross sections, dσ/dy(Pb + Pb → Pb + J/ψ + Pb). These consist of a rather non-trivial numerical evaluation of the differential NLO forward photoproduction cross sections dσ/dt(γ+Pb → J/ψ+Pb) at vanishing Mandelstam variable t based on Ref. [30], supplemented with a straightforward computation of the nuclear form factor to account for the t dependence of the cross section, as well as a nontrivial numerical evaluation of the photon fluxes from the colliding Lead nuclei based on Refs. [60,61]. In the current exploratory NLO study we adopt the simplest possible, forward-limit, approximation for the GPDs where they become just the usual PDFs. With such a "bare bones" GPD/PDF NLO framework, our goal is to test as transparently as possible, and without any additional normalization factors (which typically appear in LO studies) or modeling, how directly and efficiently the exclusive photoproduction of J/ψ mesons in Pb+Pb UPCs at the LHC actually probes the nuclear gluon distributions.
In what follows, we will first chart the scale dependence of the NLO cross sections, and compare the situation with the LO case, too. Even though the scale dependence of the NLO cross sections is known to be quite strong [30], we will show that interestingly a reasonable "optimal" scale choice can be found, with which we can, perhaps contrary to our initial expectations, simultaneously reproduce the 5.02 TeV ALICE mid-rapidity [36] and the LHCb forward-rapidity [41] data, and also the 2.76 TeV ALICE [37,39] and CMS [40] data. We will also study the corresponding NLO cross sections in photon-proton collisions, as well as their scale dependence, against the HERA and LHCb data.
We will also break down the NLO calculation into the contributions from the imaginary and real parts, as well from the gluon and quark PDFs, and show (in accordance with Ref. [30]) that the real part of the amplitude as well as the quark contributions both have a sizeable contribution and hence must not be neglected. This result indicates that, contrary to what is often claimed based on the LO results, exclusive J/ψ photoproduction in UPCs is not as direct a probe of the gluon distributions as perhaps previously thought. We will chart, by comparing the predictions obtained with the EPPS16 [44] nuclear PDFs and CT14NLO free proton PDFs [62], and nCTEQ15 [43] and nNNPDF2.0 [45] nuclear PDFs, how the gluon and quark PDFs manifest themselves in the J/ψ photoproduction UPC cross sections at different rapidities. In particular, using EPPS16, we will show that the manifestation of the nuclear effects is non-trivial and influenced especially by the relative signs of the different contributions in the amplitude. Finally, as one of the main goals of the paper, we will study how the uncertainties of the nuclear and free-proton PDFs propagate into the J/ψ photoproduction UPC cross sections.
The rest of this paper will proceed as follows: To make our study more accessible especially for the heavy-ion community and non-GPD-experts in general, we will recapitulate the theoretical NLO framework with collinear factorization and GPDs/PDFs in Sec. 2. Also the calculation of the photon fluxes and evaluation of the necessary nuclear form factors are presented there. The main results of the paper, the numerical evaluation of the co-herent exclusive J/ψ photoproduction cross sections in Pb+Pb UPCs at the LHC, their analysis and comparison with the experimental data, are presented in Sec. 3. Finally, a discussion and outlook are given in Sec. 4.

A. Differential Cross Section
In this section, we recapitulate the theoretical framework we use in our calculations for the exclusive process where A 1,2 denote the colliding nuclei and V is some vector meson (in this paper V = J/ψ). The initial-state momenta are labeled by p i and the final-state momenta by p i . Within the equivalent-photon (Weizsäcker-Williams) approximation [3,63,64], the total cross section can be expressed as where dN Ai γ (k)/dk is the centrality-integrated distribution (or flux) of photons from the nucleus A i as a function of photon energy k, and σ γ(k + )A2→V A2 and σ A1γ(k − )→V A1 are the cross sections for the photoproduction processes In the equivalent-photon approximation the photon momenta k 1,2 are considered to be collinear with colliding nuclei, and | k 1,2 | = k ± . The experimental data in Pb-Pb collisions [36,38,40,41] are differential with respect to the rapidity y of the vector meson. At fixed rapidity and transverse momentum p T of produced vector meson, the photon momentum can be expressed as where t refers to the square of the momentum transferred to the target nucleus, t = (k 1,2 − p 3 ) 2 , and M T = M 2 V + p 2 T is the transverse mass. In the typical case |t| M 2 V and p 2 T M 2 V (see e.g. Ref. [65]) so that to a very good approximation It then follows that

B. Photoproduction Cross Section
We will assume that the invariant matrix element M γA→V A for the photoproduction process can be factored into two parts, the matrix element evaluated at t = 0 and a nuclear form factor F A (t) (also called the two-gluon form factor [12]) [66], where N labels a bound nucleon and W is the c.m.s. energy of the photon-nucleon collision. It follows that the photoproduction cross section then becomes where |M γN →V N A | 2 is the square of the per-nucleon matrix element averaged (summed) over the initial-state .02 TeV. We model the form factor as the Fourier transform of the Woods-Saxon distribution [67], taking |q| = |t|. We take d = 0.546 fm [68] for the skin depth and for the nucleus radius R A we use the parametrization (see e.g. [69]), The normalization ρ 0 is fixed by requiring that F A (0) = A.
When considering the γ + p collisions we take the photoproduction cross section to be of the form [30], with [35], b GeV 2 = 4.9 + 4α P ln where W 0 = 90 GeV and α P = 0.06. This parametrization grows more slowly with W than that in Ref. [70], but is still compatible with the HERA data for exclusive J/ψ photoproduction. We have chosen the slope parameter α P to be compatible with Model 4 of [71] which fits a wider variety of elastic pp data.

C. Photoproduction Amplitude
The NLO expressions for the matrix element M γN →V N A (W, t) for photoproduction are well established in the literature [30,72] and the more recent electroproduction results [35,49,50] coincide with these in the limit of an on-shell photon. In these calculations the vector meson is considered as a composite particle of two heavy quarks in the non-relativistic approximation with zero relative velocity [73][74][75][76]. The invariant matrix element can be written as where α QED is the fine-structure constant, m Q the mass of the heavy quark, e Q the fractional charge of the heavy quark, ε V the polarization vector of the produced vector meson, ε γ the polarization vector of the incoming photon, and O 1 V is a non-relativistic QCD matrix element associated with the vector meson. Equation 13 defines the factor C which we will use later. The value of O 1 V is solved from the NLO expression for the vector-meson leptonic decay width [30,[77][78][79], where α s (µ R ) is the QCD coupling at a renormalization scale µ R . The variable ξ that appears in the Ji's parametrization of momenta [80], is the so-called skewedness parameter. In the t M 2 V limit, The function I(ξ, t) is given by where T g (x, ξ) and T q (x, ξ) are the hard-scattering coefficient functions corresponding to gluon and quark contributions [30], Here the term proportional to α s (µ R ) in T g is the purely gluonic LO contribution and the rest in T g and the whole T q constitute the NLO contributions. The exact forms of the functions f g and f q are given in Refs. [30,35,72] and we will be using specifically those of Ref. [30]. The parameter is positive and the function I(ξ, t) is understood to be evaluated in the limit → 0. Finally, F g (x, ξ, t, µ F ) is the gluon GPD and F q,S (x, ξ, t, µ F ) is the quark singlet GPD given by where µ F denotes the factorization scale. As we will consider factorization scales above the charm mass threshold, also the charm quarks are included in the above sum in conjunction with GPDs/PDFs defined in variableflavour-number schemes. As indicated in Eq. (6), we will calculate the amplitude in the approximation in which t = 0. In addition, in the current exploratory study we will approximate the GPDs by their values at ξ = 0 so that we effectively replace the GPDs with PDFs, where x ∈ [0, 1], and g(x, µ F ) and q(x, µ F ) are the gluon and quark PDFs. The differential cross section can then be written as where In our analysis we take all constants, such as the mass and the decay width of the J/ψ, from the Particle-datagroup listing [81]. The value of α s (µ R ) is taken from the LHAPDF interface [82] so that the coupling is taken consistently to be the same as the one used in defining the PDF values. The QED coupling, α QED , is evaluated throughout the work up to one loop accuracy. In our framework, following Ref. [30], we explicitly set M V = 2m Q which is an inherent assumption in our nonrelativistic approximation of the J/ψ wavefunction. In our results for |I| 2 we consistently include both the real part and the imaginary part to the results. The integrals in Eq. (21) are evaluated numerically by keeping the parameter finite but small enough so that the results are independent of its exact value. We have cross checked our numerical implementation against the method used in Ref. [35]. The factorization and renormalization scales are taken to be equal, µ = µ F = µ R , and we consider scale variation between µ ∈ [m Q , 2m Q ].

D. Photon Flux
The number of equivalent photons of energy k at a fixed transverse distance b from the center of a nucleus A with Z protons can be written as [3,13,60,61] where F is the Fourier transform of the form factor in Eq. (8) normalized to one, F (q) = F A (q)/A, and J 1 is the cylindrical modified Bessel function of the first kind. To obtain the minimum-bias flux appearing in the expression for the cross sections, e.g. in Eq. (1), we integrate over the entire impact-parameter plane multiplying N A γ (k, b) by the Glauber-type probability [83] of having no hadronic interaction, where σ NN (s) is the total (elastic + inelastic) hadronic nucleon-nucleon cross section for which we use 90 (80) mb at √ s NN = 5.02 (2.76) TeV [81], and T AA ( b) is the nuclear overlap function where T A ( b) is the nuclear thickness function, with r 2 = z 2 + b 2 and z being the longitudinal coordinate. The integrand in Eq. (23) oscillates very rapidly at large values of b, and to improve the convergence we follow Ref. [84] by making use of the flux of photons from a point-like particle. In this case one takes the nuclear density to be a delta function, ρ pl (r) = δ 3 (r), which leads to [5,85], where K 0 and K 1 are modified Bessel functions of the second kind, and The integral over the impact-parameter plane with a condition | b| > b min is also well known [86], We now rewrite Eq. (23) by adding and subtracting the flux of photons from a point-like particle, By taking b min = 30 fm or higher, the last term will be negligible. Differences between this result and the pointlike approximation have been studied e.g. in Refs. [3,84].

A. Absolute magnitude and scale sensitivity of cross sections
First, we chart the uncertainty arising from the choice of the factorization/renormalization scale in the exclusive rapidity-differential J/ψ photoproduction cross sections NLO with EPPS16 s NN = 2.76 TeV The scale-choice uncertainty-envelope of the rapidity-differential exclusive J/ψ photoproduction cross section in ultraperipheral Pb+Pb collisions at √ sNN = 5.02 TeV, as a function of the J/ψ rapidity y, calculated to NLO pQCD with the EPPS16 nPDFs [44] and compared with the experimental data from Refs. [38] (ALICE Forw), [36] (ALICE Cent) and [41] (LHCb Forw). The experimental data points are mirrored w.r.t. y = 0, and their errorbars are obtained by adding the statistical and systematic errors in quadrature. The solid (red) curve shows the NLO result with our "optimal" scale explained in the text. Lower panel: The same but at √ sNN = 2.76 TeV and with experimental data from Refs. [39] (ALICE Forw), [37] (ALICE Cent) and [40] (CMS Cent). For the errorbars of the data, all given errors are added in quadrature.
in ultraperipheral Pb+Pb collisions. Figure 1 shows the uncertainty envelopes that result from varying the scale (upper panel) and 2.76 TeV (lower panel), using the central set of the EPPS16 nPDFs [44]. For comparison, the figure also shows the experimental LHC data measured at these energies at forward rapidities by ALICE [38,39], LHCb [41] and CMS [40], and at central rapidities by AL-ICE [36,37]. The solid (red) lines in the middle-parts of the envelopes show the results with µ = 0.76M J/ψ = 2.37 GeV, a scale we have iteratively obtained by requiring a rough simultaneous fit to the data at both collision energies. In what follows, we call this "optimal" scale, emphasizing however that its precise number bears no special significance but it depends e.g. on the assumed the GPD modeling details and nPDFs in general.
On the one hand, as expected based on Ref. [30], we observe that the scale uncertainty remains quite large also here in the nuclear case. On the other hand then, it is interesting and quite encouraging that already with our current "bare bones" GPD/PDF framework the NLO cross sections with entirely feasible scale choices µ =  [62] PDFs and compared against the experimental HERA data from H1 [26] and ZEUS [27], and LHC data from LHCb [28,29]. The solid (red) line corresponds to the "optimal" scale explained in the text.
O(M J/ψ ) not only are of the correct order of magnitude but actually some scale-choices can be found with which we can rather well reproduce the data at all rapidities and both collision energies. Earlier, especially with (ad hoc normalized) LO cross sections and the forward ALICE data at 5.02 TeV, this seemed not to be the case [36].
Second, as a further check of our UPC results from the "bare bones" GPD/PDF framework, we study in Fig. 2 to what extent we can reproduce the exclusive J/ψ photoproduction cross sections measured in ep collisions at HERA and in pp collisions at the LHC. 1 The NLO cross sections here are, for consistency, computed with the CT14NLO PDFs [62] which is the free-proton PDF set that the EPPS16 nPDFs are based on. The envelope shows again the uncertainty arising from varying the scale µ between M J/ψ /2 and M J/ψ . The HERA data in the figure are from H1 [26] and ZEUS [27], and the LHC data from LHCb [28,29]. The solid (red) line in the middle of the envelope is again the NLO cross section computed with our "optimal" scale which reproduced the nuclear data. As expected based on Ref. [30] and other previous NLO studies of this process [31][32][33][34][35], the scale dependence is indeed large, and especially towards larger values of the photon-proton c.m.s. energy W the data easily fall within the envelope. From the point of view of the nuclear UPCs the most relevant c.m.s.-energy region here is W = 10 . . . 700 GeV (see the second x-axis in Fig. 5 ahead). Interestingly, our framework with the 1 The photoproduction cross sections are extracted from the LHC pp data through rather minimal modeling [28,29].
"optimal" scale leads to a rather reasonable overall agreement with the HERA/LHC ep/pp data as well, except perhaps for the very lowest W points. As suggested by earlier work [30][31][32][33][34][35], there is room for GPD modeling testable against the ep/pp data but given the large scaleand PDF-uncertainties (discussed in Fig. 13 ahead), and also the exploratory nature of the present NLO study for UPCs of nuclei, we leave this as future improvement. With Fig. 2, it is also worth emphasizing that in the previous LO UPC studies one has typically normalized the LO cross sections to the HERA/LHC ep/pp data and carried the obtained normalization factor then over to the UPC study, while in our current NLO study there are no ad hoc normalization factors.
Third, we investigate the stability of the rapiditydifferential J/ψ photoproduction cross sections in Pb+Pb UPCs, i.e. the changes in the magnitude and shape, and in the scale-dependence of the cross sections, when moving from LO to NLO in pQCD. These questions are answered by Fig. 3, where we show the rapiditydifferential cross sections computed with various fixed scales µ between M J/ψ /2 and M J/ψ in the LO and NLO cases (upper and lower panels, correspondingly). To be exact, the LO here refers to the purely gluonic Bornterm contribution which enters the full NLO result. For the computation, we again use the EPPS16 nPDFs. We observe that the overall effect of the NLO terms is to reduce the LO cross sections rather significantly, at the "optimal" scale by a factor of 2.3 at mid-rapidity, and by a factor of 3.3 at y = ±4. We also see that the stud- tations based on Ref. [30] also now in the nuclear UPC case, that at the low scales of µ = O(M J/ψ ) the NLO contributions do not stabilize the results, yet, but bring the cross sections nevertheless into the right direction. Interestingly, as seen in Fig. 3, also the whole shape of both the LO and the NLO results is quite sensitive to the scale µ, and again perhaps even more so at NLO, in this scale range. In the LO case, the strong scale dependence can be traced back mainly to the rapidly changing gluon distributions, while in the NLO terms the scale µ resides both in the pQCD matrix elements and in the PDFs. In particular, as we will soon see, in the NLO cross sections the rapidly evolving small-x quark PDFs start to play a surprisingly important, and at mid-rapidities even a dominant, role.
To analyse the scale dependence of our LO and full NLO results and their interrelation further, we plot in Fig. 4 the computed rapidity-differential cross sections at fixed rapidities y = 0 and at y = ±4 as a function of the scale µ. As we see, at y = 0 the scale dependence at low scales is stronger in the NLO than in the LO results but towards higher scales it actually becomes weaker. At y = ±4 we see the scale dependence being stronger in NLO at all scales studied. Thus, whether the scale dependence is improved (tamed) when going from LO to NLO depends on the rapidity y and potentially also the scale-choice region. Another interesting observation is that our "optimal" scale µ = 2.37 GeV is right in the region where the scale dependence at y = 0 turns from stronger to weaker relative to LO, i.e. where the LO and NLO results are closest to each other. At y = ±4, however, we do not find a similar taming effect to take place. This figure also shows how the NLO/LO ratio ("K-factor") is not a constant as a function of the scale, and certainly not a constant as a function of the J/ψ rapidity.
B. Complex structure of the cross section Next, we discuss the very interesting consequences of the complex structure of the rapidity-differential J/ψ photoproduction cross sections in 5.02 TeV Pb+Pb UPCs. First, in Fig. 5 we study the k ± contributions in Eq. (4) to the rapidity-differential J/ψ photoproduction NLO cross section in 5.02 TeV Pb+Pb UPCs, computed with EPPS16 at our "optimal" scale. The photon-proton c.m.s. energy corresponding to the photon energies k ± in Eq. (4) are denoted by W ± in what follows. As indicated by the second x-axis at the top of Fig. 5, W + (W − ) increases to the right (left). As we saw in Fig. 2 above, the photoproduction cross section in the k ± terms of Eq. (4) increases as a function of W ± , correspondingly. The photon flux, however, decreases rapidly as a function of the energy W ± (see e.g. Fig. 3 of [61]), causing the nonmonotonous behaviour of the two symmetric contributions as is seen in Fig. 5. Looking at the W + curve (dashed, red) we see that first at backward-most rapidities the photon flux is high enough to produce a noticeable cross section in spite of the smallness of the photoproduction cross section there. Also the t integral of the squared form factor of Eq. (8) reaches non-negligible values by y ∼ −4, which also contributes to the initial rise of the cross section at backward-most rapidities. Then in the "shoulder" region the decrease of the photon flux wins over the increase of the photoproduction cross section, causing the small dip seen in the figure. Approaching then mid-rapidities, the increase of the photoproduction cross section now wins over the decrease of the photon flux, until eventually towards forward-most rapidities the photon flux decrease again dominates and the resulting cross section dies out. For the W − component (dotted green curve), the behaviour is a mirror image of this, and the final result (solid blue curve) is a combination of the W ± contributions as seen in the figure.
Second, we quantify the contributions from the imaginary and real parts of the amplitude. The decomposition of the full result (∝ |M| 2 ) into the contributions from the real part (∝ |Re(M)| 2 ) and the imaginary part (∝ |Im(M)| 2 ) for both the LO and NLO cross sections is shown in Fig. 6. These results are again obtained with the EPPS16 nPDFs and fixing µ to our "optimal" scale. The LO here again refers to the Born term contributions entering the full NLO result. As the upper panel shows, in the LO case where only gluons contribute, we confirm -at least for gluon PDFs of a modest small-x rise, such as those in EPPS16/CT14NLO -the general claim that the contribution from the imaginary part of the amplitude clearly dominates at all rapidities. However, as the lower panel shows, the situation changes rather dramatically for the NLO cross sections: At mid-rapidity the contribution from the real part of the amplitude is about a quarter, which clearly is not anymore negligible. Towards forward/backward rapidities the real-part contributions become even more important and, as seen in the figure, there is a region at large/small rapidities where they dominate over the imaginary-part contributions. These findings are also consistent with those of  (4) to the NLO exclusive rapidity-differential J/ψ photoproduction cross section in 5.02 TeV Pb+Pb UPCs as a function of the J/ψ rapidity y, computed using EPPS16 nPDFs and with our "optimal" scale. The second x axis on the top shows the values of W + corresponding to each y.
Ref. [30], see Fig. 17 there. The message from Fig. 6 is clear: both the imaginary and real parts of the amplitude must be accounted for in the calculation of these cross sections. Third, in Fig. 7, we investigate the breakdown of the computed J/ψ photoproduction NLO cross section in 5.02 TeV Pb+Pb UPCs into the quark and gluon contributions, using EPPS16 and our "optimal" scale. The solid (blue) curve labelled "Full |M| 2 " is the full NLO cross section of Fig. 6, while the dashed red (dotted green) curve labelled "Only Gluons" ("Only Quarks") is obtained by setting the quark (gluon) distributions to zero. The dashed-dotted curve labelled "Interference" corresponds to the remaining contribution from the cross section pieces that contain both quarks and gluons. As shown by Fig. 7, at mid-rapidity the quarks-only contribution dominates over the gluons-only by a factor of four (!), and the quark-gluon term over the gluons-only by a factor of three. Towards forward/backward-most rapidities the gluons-only contributions become the dominant ones, and we can see that the gluon-quark term also changes its sign when going from mid-rapidity to forward/backward rapidities.
Recalling the original attraction of the exclusive J/ψ photoproduction in electron-proton collisions and in nuclear UPCs as an exceptionally efficient probe of small-x gluon distributions, the results in Fig. 7 appear at first sight somewhat surprising. Especially the quark dominance at mid-rapidity seems to be in a direct contradiction with the original LO-based gluon-probe suggestion, and in fact also with our expectation that small-x gluons should after all dominate also the NLO contributions! A better understanding of this clearly calls for a more so that the squared amplitude entering the cross section becomes The gluons-only contribution in Fig. 7 comes from the term (33) and the quarks-only contribution from while the gluon-quark interference contribution corresponds to the third term on the r.h.s. of Eq. (32). Figure 8 shows the above real and imaginary parts of the amplitude, multiplied with the factor ξ/C (see Eq. 13), as a function of the rapidity corresponding to W + . ) changes its sign from plus to minus when approaching backward rapidities, which causes the sign change of the quark-gluon mixing term in Fig. 7. Let us look at the following three example-rapidities: • At y = 0, where the W ± components contribute equally (see Fig. 5), the cancellation of the gluon terms is coincidentally (that is, with these PDFs) almost perfect in the imaginary part, so that which makes the imaginary part of the quark amplitude dominate the cross section in Fig. 7. In the quark-gluon mixing term then the product of the imaginary parts dominates over the product of the real parts, and due to the large Im(M NLO Q ) the quark-gluon contribution dominates over the gluons-only term.
• At y ≈ −3, Fig. 5 indicates that the W ± contributions are equally important, so that Fig. 8 should be read both at y ≈ −3 and y ≈ +3. The squared amplitude is larger for y = 3 but the rapid decrease of the W − - component's photon flux and nuclear form factor towards negative rapidities now suppresses the W − component so that it becomes of the same magnitude as the W + component whose squared amplitude is smaller but photon flux correspondingly larger. As seen in Figs. 6 and 7, as a result of these competing effects the real and imaginary parts of the amplitude, as well as quarks and gluons, then contribute equally to the rapidity-differential cross section at y ≈ −3.
• At y ≈ −4, where the cross section is dominated by the W + component as seen in Fig. 5, the LO and NLO gluon terms cancel to a much smaller degree both in the real and imaginary parts, and the hierarchy becomes causing the gluons-only terms to dominate over the quarks-only by a factor of four. In this case, the sizable quark-gluon mixing term is deeply negative because of the large negative term Im(M LO G ) + Im(M NLO G ). It is again the negative sign of this term that in the full amplitude causes [Re(M)] 2 [Im(M)] 2 , seen in Fig. 8 and in the lower panel of Fig. 6 at y = −4... − 3.
As shown by Figs. 5-8, the full NLO cross section thus has a very detailed complex structure with interplays between the photoproduction cross section, the photon flux and the nuclear form factor, between the W ± components, and especially between the various contributions from the real and imaginary parts of the amplitude.
The key to understand the obtained rapidity-differential cross sections is the degree of cancellation of the LO and NLO gluon contributions of opposite signs. We have also checked that the situation is qualitatively the same for the 2.76 TeV collision energy, and that the real part contributions become slightly more important for all values of y than for the 5.02 TeV case. We have also checked that in the case of no nuclear effects, the situation remains qualitatively the same.

C. Nuclear effects and PDF uncertainties in the cross section
Next, we analyse how the nuclear modifications of the PDFs as well as the uncertainties of the nuclear and free-proton PDFs propagate into the exclusive rapiditydifferential J/ψ photoproduction cross sections. Figure 9 compares the rapidity-differential cross sections at 5.02 TeV obtained at our "optimal" scale with the EPPS16 nPDFs (solid orange curve), and the one obtained with the CT14NLO free-proton PDFs (dashed blue) which are the baseline for EPPS16. As seen in the figure, at midrapidity, where the W ± terms contribute equally, the cross sections show a reduction of a factor of 0.76 from CT14NLO to EPPS16. Towards backward/forward rapidites, i.e. in the regions where the W ± terms contribute significantly and probe the nuclear effects in different x-regions, the net nuclear effects are slighly increasing. Finally at the backwardmost (forwardmost) rapidities, where the single W + (W − ) contribution dominates and one enters the antishadowing region, the nuclear effects essentially die out.
The general behaviour and magnitude of the nuclear effects here can be understood as follows: • First, we recall from Figs. 7 and 8 that it is the imaginary part of the quark amplitude that dominates the cross section at y = 0. Recalling that ξ(y) = ζ(y)/(2 − ζ(y)), where ζ(y) = M 2 J/ψ /W 2 and W 2 = M J/ψ e y √ s NN , we have ξ(y = 0) ≈ 3 · 10 −4 . This is deep in the shadowing region of nPDFs, and in EPPS16 at this x and our "optimal" scale the average nuclear sea-quark (gluon) modification is about 0.68 (0.74). The fact that there seems to be a weaker than quadratic dependence on the PDF's nuclear modification factor follows to our understanding from two reasons: First, in the NLO amplitudes one integrates the parton distributions over x from zero to one: At x ξ shadowing is about a constant factor (in EPPS16) while at x ξ shadowing diminishes, so that the net effect of the x-integration is a taming of the nuclear effect from that at x = ξ. Second, and most importantly, as discussed in detail below, it is again the surprisingly complicated interplay of the different parts of the amplitude and in particular the mutual cancellation of the LO and NLO gluon amplitudes that causes the quark-gluon mixing term to actually cancel some of the nuclear effects.
• Towards backward rapidities there are competing nuclear effects as W + decreases, the probed values of ξ increase and the nuclear modifications thereby decrease, and as simultaneously W − increases, the probed values of ξ decrease and the nuclear modifications thereby increase (and towards forward rapidities conversely). And, as seen in Fig. 7 also quarks and gluons compete over the dominance of cross section, the quark dominance turning into a gluon one towards backward/forward rapidities.
• At the backward rapidity y = −4 then, we recall the W + contribution (Fig. 5) and the NLO real part of the full amplitude (Figs. 6 and 8) to dominate the cross section, and from Fig. 7 we again see that both quarks and gluons contribute here. Now ξ(y = −4) ≈ 1.7 · 10 −2 and the EPPS16 gluon (quark) modification is a factor of 0.88 (0.86) while the net nuclear effect is about a factor of 0.68, i.e. surprisingly large. In this region the integration over x does not tame the nuclear effects to the same degree as at small values of x, and in particular the large and negative quark-gluon mixing term drives the efficiency of nuclear effects up here.
Given the complex intertwined structure of the cross section, it is also useful to analyse what happens if we start from the EPPS16 result and turn separately off the nuclear effects from gluons and quarks, one at the time.
• Turning first off the nuclear effects (suppression) in the gluon PDFs results in the dashed-dotted (green) curve labelled "Gluons with CT14NLO" in Fig. 9, which shows a reduction in the cross section relative to the EPPS16 result (solid orange curve) at mid-rapidity. This seems again quite counter-intuitive, as we would naively expect a removal of suppression to cause an increase instead. Such a behaviour can, however, be again understood by studying the real and imaginary parts of the amplitude: ) in the quark-gluon mixing term, causing the suppression that we see in Fig. 9 at mid-rapidity.
• Then, turning off the nuclear effects in the quark distributions, but leaving them on in the gluon contribution results in the dotted black curve which lies rather close to the pure CT14NLO case of no nuclear PDF effects at all. This time this is an obvious result, as at mid-rapidity the quark part Im(M NLO Q ) dominates the cross section and removing the suppression in the PDFs just increases the cross section as expected. Figure 9 thus underlines the quark dominance demonstrated earlier in Fig. 7.
Because of the rather counter-intuitive results above, and since there is the integration over x from zero to one in the NLO amplitude, we would like to confirm that NLO exclusive photoproduction of J/ψ in Pb+Pb UPCs at the LHC indeed probes the small-x shadowing region (x 0.03...0.04 in EPPS16), and not the antishadowing region (0.03...0.04 x 0.3 in EPPS16) in the nPDFs. If the process indeed probes the quark and gluon distributions at x = O(ξ), and ξ(y = 0) ≈ 3 · 10 −4 , then the biggest effect to the final result (relative to the CT14NLO result above) should be attained by turning on only the nuclear corrections in the shadowing region. We have checked that this is indeed the case: Running the code with ad hoc modified nPDFs that coincide with EPPS16 in the shadowing region and with CT14NLO elsewhere, the results are essentially (within 6%) the same as the EPPS16 results.
Next, we investigate how sensitive the studied cross sections are to the choice of the nPDFs. Figure 10 shows the rapidity-differential cross sections obtained with the central sets of the EPPS16 (solid orange curve), nCTEQ15 (dashed green) and nNNPDF2.0 (dotted blue) nPDFs. The nCTEQ15 set gives essentially the same result as EPPS16 but there seems to be a huge difference to the nNNPDF2.0 set. The shape of the nNNPDF2.0 result is very different from EPPS16/nCTEQ15, and the magnitude at forward and backward rapidities is off by about a factor of 15. We have traced the very fast growth of the cross section down to the rapidly growing real part of the LO gluon amplitude, which includes again the integration over x from 0 to 1 where the small-x gluons (in the ERBL region x ξ but near x ∼ ξ) start to play a significant role with nNNPDF2.0. The real part of the LO gluon amplitude is not as well numerically cancelling against the real part of the NLO gluon amplitude with nNNPDF2.0 as with EPPS16/nCTEQ15, which in turn makes the forward/backward-y cross section again more sensitive to the small-x gluon distributions, and this is what we see in Fig. 10.
We plot in Fig. 11 the gluon distributions xg(x, µ) and the quark singlet distributions F q,S = q [q(x, µ) + q(x, µ)] from EPPS16, nCTEQ15 and nNNPDF2.0 nPDFs as they enter our computation at the "optimal" scale. The figure confirms the similarity of the EPPS16 and nCTEQ15 PDFs and shows that the nNNPDF2.0 quarks differ systematically from these at x 10 −5 and the gluons at x 10 −4 . In Fig. 10, the increased small-x gluons of nNNPDF2.0 make the W − component of the cross section to be the dominant one at y = −3. For the W − contribution ξ(y = −3) = O(10 −5 ), and at these values of x, Fig. 11 indicates already a factor of three difference between the nNNPDF2.0 and EPPS16/nCTEQ15 gluons. The square of this difference then explains the order of magnitude of the difference between the nNNPDF and EPPS16/nCTEQ15 results seen in Fig. 10.
Next, we investigate the PDF uncertainties in the computed rapidity-differential cross sections and compare them with the existing data. We propagate the PDF/nPDF uncertainties to the computed cross sections using the asymmetric form [44] where S ± i labels the error sets for the given PDF. We plot the error sets of EPPS16+CT14NLO in Fig. 12 for the gluon distributions xg(x, µ), and for the quark singlet distributions F q,S , again at our "optimal scale". As the figure shows, one CT14-related error set, "Set93", of the EPPS16 implementation in LHAPDF [82] (error set 53 in CT14NLO), stands clearly out at smallest values of x, and even more strongly than the nNNPF2.0 PDFs did in Fig. 11, while the rest of the EPPS16-related and CT14related error sets show only rather moderate variations w.r.t. the central sets. Similarly to the case with the nNNPDF2.0 nPDFs above, the rapid growth of the smallx gluon distributions in this error set induces again a rapid growth of the real part of the LO gluon amplitude, and hence the cross sections. Figure 13 shows the uncertainties that are induced to the rapidity-differential exclusive J/ψ photoproduction cross sections in 5.02 TeV (upper panel) and 2.76 TeV Forward limit GPDs = F = 2.37 GeV F q, S multiplied with 10 2 11. The nPDF gluon distributions xg(x, µ) and the quark singlet distributions F q,S = q [q(x, µ) +q(x, µ)] as given by EPPS16 (solid lines), nCTEQ15 (dashed) and nNNPDF2.0 (dotted) nPDFs at the "optimal" scale.  (lower panel) Pb+Pb UPCs by the PDF/nPDF uncertainties. The uncertainties arising from the EPPS16 nuclear effects alone are shown by the dark (blue) bands, while the full EPPS16+CT14NLO error bands (green) contain uncertainties from both the nuclear effects and the free-proton baseline PDFs. The results with EPPS16 and CT14NLO central sets are shown by the solid (blue) curves. As expected based on Fig. 12, "Set93" above entirely dictates the green error bands. The EPPS16+CT14NLO full uncertainty band at mid-rapidity (not shown in the figure) goes up to some 150 (37) mb and at y ≈ ±2.2 as high as 1500 (170) mb for the 5.02 (2.76) TeV collision energy. We also have checked that without "Set93" the CT14NLO uncertainties become of the same order and slightly smaller than those for EPPS16. For comparison with the EPPS16 results, we also plot the uncertainty bands (hatched) arising from the nCTEQ15 error sets. These now account for the uncertainties in the nuclear effects only, and not in the free-proton PDFs. The central-set results with nCTEQ15 are shown by the dashed (red) line. We should also emphasize that the nCTEQ15 results here have been obtained at our "optimal" scale, without further tuning of the scale. As we have already seen, the EPPS16 results produce a relatively good fit to the experimental Run1 and Run2 data at our "optimal" scale, and as seen in Fig. 13, so do the nCTEQ15 ones, too. The uncertainties arising from the nuclear effects in EPPS16 and nCTEQ15 are of the same order of magnitude mutually, and typically somewhat larger than the errorbars of the data. As the figure indicates, one must not forget the free-proton PDF uncertainties when considering absolute cross sections. Finally, regarding the tension between the ALICE and LHCb data in the forward/backward direction, we can see that at least at our "optimal" scale both the EPPS16 and nCTEQ15 results (but obviously not the nNNPDF2.0) seem to reproduce the LHCb data points better but that both data sets can still be accommodated within the larger EPPS16 uncertainties.

IV. SUMMARY
We have presented the very first implementation of exclusive rapidity-differential J/ψ photoproduction cross sections in ultraperipheral nucleus-nucleus collisions in the framework of collinear factorization and NLO per-turbative QCD. We have developed our numerical code for the ultraperipheral nuclear collisions based on the analytical NLO results of Ref. [30], utilizing the experience obtained also in [35,48], and following earlier literature in accounting for the photon fluxes of the colliding nuclei [3,5,13,60,61,[84][85][86] and for the t-dependence of the cross section with a standard nuclear form factor. In this exploratory NLO study for the UPCs, we approximate the GPDs involved in the process with their forwardlimit nuclear PDFs. Our default choice for the nPDFs and their error sets is EPPS16 [44] but we also study the nPDF sensitivity of our results by using nCTEQ15 [43] and nNNPDF2.0 [45].
We have shown that, as expected based on Ref. [30], the computed rapidity-differential NLO cross sections of J/ψ photoproduction in 5.02 and 2.76 TeV Pb+Pb UPCs at the LHC, as well as the corresponding photoproduction cross sections in ep collisions at HERA, are both in their magnitude and in their shape quite sensitive to the scale choice. As the scale-sensitivity is much larger than the error bars of the experimental data at the LHC, it makes it difficult to make solid NLO predictions of the corresponding J/ψ cross section for UPCs at other energies. Quite encouragingly, however, we have found that a scale-choice µ ≈ 0.76M J/ψ , which lies in the physically reasonable range µ = O(M J/ψ ), can actually be determined, with which we can well reproduce the ALICE [36][37][38][39], LHCb [41] and CMS [40] UPC data at these energies. We have also tested that the same scale choice, called here "optimal" scale, works well also with the nCTEQ15 nPDFs. Interestingly, in studying the scale-sensitivity at a fixed value of y = 0, we noticed that towards the upper end of the scales studied here the scale sensitivity of the full NLO result becomes actually weaker than that of the LO result, but towards the lower end of scales stronger than in LO. Also interestingly, at midrapidity the "optimal" scale becomes fixed right in the scale-region where the NLO contributions are the smallest relative to LO. In the future, it will be interesting to see whether this "minimal-sensitivity" feature remains there also after further modeling of the GPDs.
We have made an effort to analyse in sufficient detail the surprisingly complex structure of the exclusive rapidity-differential J/ψ photoproduction NLO cross sections in Pb+Pb UPCs at 5.02 and 2.76 TeV. In particular, we have shown how the computed NLO cross sections form under various competing and intertwining effects: There are competing contributions from the photon-nucleon c.m.s. energy W ± components, from the real and imaginary parts of the full amplitude, from the quark and gluon GPD/PDF contributions which also mix in a non-trivial way in the squared amplitude, and most importantly of all, from the gluonic LO and NLO amplitudes which come with opposite signs and cancel each other to a degree that nontrivially depends on the W ± . All these competing contributions need to be taken into account in the full NLO study, as is done in the current paper.
The main result of our NLO study with the EPPS16 nPDFs, similar to the findings in Ref. [30] but now for UPCs, is that due to the cancelling LO and NLO gluon amplitudes it is predominantly the small-x quark GPDs/PDFs that exclusive J/ψ photoproduction is probing in UPCs at midrapidity, and not the gluon distributions as has been traditionally suggested before based on LO. This is an important result not addressed to our knowledge in the UPC context before. We have also checked that this result is robust against the scale variation studied here. We have also shown that towards the forward/backward rapidities the gluon dominance is eventually recovered but because of the folding with the photon flux (which kills one of the W ± contributions) the nuclear gluon GPDs/PDFs become probed at larger values of x (where shadowing effects become smaller) than at midrapidity. Thus, our conclusion is -at least in our current "bare bones" GPD/PDF framework and with the EPPS16 and nCTEQ15 nPDFs -that the exclusive rapidity-differential J/ψ photoproduction cross sections at the LHC are not as a direct and efficient probe of the small-x nuclear gluon PDFs as thought before, but that they are primarily probed (at midrapidity at least) through the DGLAP evolution of the quark GPDs/PDFs. Another important observation is that at midrapidity the dependence of the computed NLO cross sections on the nuclear effects in PDFs is not as quadratic as thought in the LO gluon context before. The taming of the net nuclear effects follows partly from the x-integration in the NLO amplitude but predominantly from the behaviour of the interference term in the squared amplitude which mixes the quark and gluon contributions in a non-trivial way.
We have also investigated the dependence of our results on the uncertainties of the PDFs. The nCTEQ15 centralset results are essentially the same as those with the central set of EPPS16. At midrapidity, where the quark contributions dominate, these two sets show very similar error bands when the uncertainties of the nuclear effects in the PDFs are propagated into the NLO cross sections. Towards forward/backward rapidities where gluons dominate, the EPPS16 uncertainties become slightly larger, which follows from the more realistic (due to having more freedom in the gluon PDF shape there) estimates of the gluon nPDF uncertainties than in nCTEQ15. In any case, in the current "bare bones" GPD/PDF framework, we observe that both the forward ALICE [38] and LHCb [41] data can be accommodated within the nuclear PDF error bands, while the results with the central sets of EPPS16 and nCTEQ15 agree better with the LHCb data.
Finally, we have observed that if there is a very rapid rise in the small-x gluon distributions, such as in the nNNPDF2.0 central set and the error set 53 in the CT14NLO free-proton PDFs [62] (93 in EPPS16 at LHAPDF), then the smallest-x contribution to the real part of the gluon LO amplitude starts to dominate the cross sections. Concretely, in our results when the EPPS16 nuclear errors and the CT14NLO errors are ap-propriately combined, the CT14NLO error set 53 (93 in EPPS16 at LHAPDF) dictates the upper boundaries of the very large uncertainty band on our central result. In our "bare bones" GPD/PDF framework, such a growth seems to be ruled out by the UPC data considered here. However, before we can make any further conclusions on this point, uncertainties arising from the modeling of GPDs should be quantified.
The current paper is meant as a baseline for systematic further studies of exclusive photoproduction of vector mesons in ultraperipheral nucleus-nucleus collisions, in collinear factorization and NLO pQCD. An obvious next task is to repeat the NLO study for the photoproduction cross sections of Υ mesons, to investigate in particular how much the scale dependence changes and check exactly what happens with all the intertwined effects at the higher scales µ = O(M Υ ). On the basis of Ref. [30], we would expect to see a reduced scale-sensitivity and a stronger dependence on the gluon PDFs also in the UPC case.
There are also several ways the current framework could and should be improved. Our strategy for the current exploratory study is that as the scale-and PDFrelated uncertainties are so large, we may leave the GPD modeling (such as in Ref. [55]) as a future challenge. Next, given the studied "bare bones" GPD/PDF baseline, it will be interesting to study how the nPDF uncer-tainties propagate to the GPDs and via them to the NLO cross sections. As far as we can see, based e.g. on Refs. [53,87] the skewedness corrections to the GPD quark distributions in the DGLAP region can be expected to be larger for quarks than for gluons, which would further strengthen our conclusion of the quark dominance at mid-rapidity. Towards forward/backward rapidities, the gluon dominance would then correspondingly kick in more slowly. Particularly interesting here would be to study the role of the nuclear effects in the ERBL region, where the PDFs are known not to be an optimal approximation but which in the current study turned out to be important essentially only with PDF sets that have rapidly growing small-x distributions. Future improvements would also include non-relativistic QCD corrections into the vector meson wavefunction [88][89][90].