Predictions for exclusive Υ photoproduction in ultraperipheral Pb + Pb collisions at the LHC at next-to-leading order in perturbative QCD

Wepresentpredictionsfortherapidity-differential cross sections of exclusive Υ photoproduction in ultrape-ripheral collisions (UPCs) of lead ions at the Large Hadron Collider (LHC). We work in the framework of collinear factorization at next-to-leading order (NLO) in perturbative QCD, modeling the generalized parton distributions (GPDs) through the Shuvaev transform of nuclear parton distribution functions (nPDFs). Our direct NLO predictions depend sig-niﬁcantly on nPDF uncertainties and on the choices of the factorization and renormalization scales, but to a much lesser degree on GPD modeling. To tame the scale dependence and to account for the fact that the NLO calculations generally underpredict the photoproduction measurements on protons, we also present alternative, data-driven predictions. In this approach the underlying photoproduction cross sections on lead are found by combining their nuclear modiﬁcations calculated at NLO with the measured photoproduction cross sections on protons. The data-driven strategy reduces the uncertainties associated with the scale choices, and essentially eliminates the effects of GPD modeling thereby leaving the cross sections sensitive mainly to the input nPDFs. Our estimates indicate that the process is measurable in Pb + Pb collisions at the LHC.


Introduction
The exclusive production of heavy vector mesons V in ultraperipheral collisions (UPCs) of heavy nuclei, A 1 + A 2 → A 1 + V + A 2 , has for a long time captured the interest of both the theoretical and experimental high-energy physics communities.It allows one to study not only the perturbative aspects of Quantum Chromodynamics (QCD) but to also probe the non-perturbative structure of nuclei [1][2][3][4].These ultraperipheral events are largely initiated by electromagnetic interactions as the short-range hadronic interactions are strongly suppressed by the exclusivity of the final state.The vector meson production then effectively proceeds through an interaction of a quasi-real photon from one nucleus with the other nucleus such that the colliding nuclei remain intact and the exclusivity of the vector meson production is maintained via a net-colourless production mechanism.At leading order (LO) in perturbative QCD (pQCD) [5], the production is mediated through a two-gluon exchange, while at nextto-leading order (NLO), there is also a quark-pair initiated contribution [6].The exchanged partons carry different longitudinal momentum fractions depending on an additional off-forward skewness parameter, ξ .This results in a factorization [7] of the scattering amplitude into the perturbatively calculable hard-scattering part and non-perturbative generalized parton distribution functions (GPDs) [8][9][10].
The first UPC measurements of exclusive J/ψ mesons came from the PHENIX collaboration at the Relativistic Heavy Ion Collider (RHIC) in Au + Au collisions at the nucleon-nucleon centre-of-mass system (c.m.s.) energy of √ s N N = 200 GeV [11].Subsequently, the ALICE, CMS, and LHCb collaborations at the Large Hadron Collider (LHC) have measured the same process in heavier Pb + Pb UPCs at √ s N N = 2.76 and 5.02 TeV in a wide range of the J/ψ rapidities from y = 0 up to |y| ∼ 4.5 [12][13][14][15][16][17][18].These data -not forgetting the multitude of statistics anticipated in the heavy-ion programme of the High Luminosity LHC [19] -provide ample grounds for understanding the perturbative structure of QCD and the nuclear shadowing phenomenon encoded e.g. in nuclear parton distribution functions (nPDFs) [20][21][22][23], down to momentum fractions of x ∼ (M V / √ s N N ) exp(−|y|) ∼ 10 −5 at resolution scales , where M V is the mass of the vector meson.In our previous works [24,25], we studied the exclusive photoproduction of J/ψ mesons in Pb + Pb and O + O collisions to next-to-leading order (NLO) in pQCD.By approximating the GPDs with PDFs, we demonstrated the complicated interplay of the quark and gluon contributions at NLO over the entire LHC acceptance in rapidity, and showed that our theoretical predictions agree with the experimental data for this process [12][13][14][15][16][17][18] within the large theoretical uncertainties associated with the choice of the factorization/renormalization scales and nPDFs.Valuable and complementary information on nPDFs at small momentum fractions x, in particular, on the scale dependence of nuclear shadowing, can be obtained by studying exclusive photoproduction of heavier vector mesons such as Υ mesons consisting of a bottom quark and its antiquark.To date, while there have been measurements of the exclusive photoproduction of Υ in e + p collisions at Hadron Electron Ring Accelerator (HERA) [26][27][28], as well as in p + p [29] and p + Pb collisions at the LHC [30], there has been no reported measurement of the exclusive production of Υ mesons in heavy-ion collisions.
In the work presented here, we make predictions for the rapidity-differential cross sections of this process at √ s N N = 5.02 TeV in Pb+Pb collisions, extending our previous framework to incorporate a more careful GPD modeling by relating the nPDF to nuclear GPDs through the so-called Shuvaev integral transform [31][32][33].Despite the larger interaction scale in comparison to the J/ψ case, the theoretical uncertainties in the case of Υ production are still sizable and -as was already noticed in the pioneering work of Ref. [6] and as we confirm in this work as well -natural choices of the factorization/renormalization scales μ 2 ∼ O(M 2 Υ ) do not lead to a particularly good description of the HERA data.This would then cast doubts also on our direct NLO predictions in Pb + Pb.As a workaround, we will adopt an alternative method in which we anchor our predictions for the underlying γ +Pb → Υ +Pb cross sections on the HERA data on the γ + p → Υ + p process by using the NLO calculations only for the ratios of cross sections between these two processes.We call this method the data-driven approach.We also analyze the nuclear modifications of the γ + Pb → Υ + Pb cross sections due to nuclear effects in PDFs and show that for ξ < 10 −3 , they coincide very closely with the gluon nuclear modification factor squared.
The rest of the paper is organised as follows.In Sect.2.1, we summarize our theoretical framework for the exclusive photoproduction of Υ in ultraperipheral Pb + Pb collisions within NLO pQCD, and then discuss the modeling of GPDs in Sect.2.2.The ingredients of our data-driven approach are explained in Sect.3. In Sect.4, we then present our results for the cross sections and their nuclear modifications, discussing also how our calculations build up from various components.Finally, we draw our conclusions in Sect.5, outlining also future directions.

Theoretical framework
2.1 Exclusive Υ photoproduction in Pb + Pb UPCs at NLO pQCD Within the equivalent-photon approximation [1,2], the rapiditydifferential cross section for the process Pb + Pb → Pb + Υ + Pb can be written as where kd N Pb γ (k)/dk is the Weizsäcker-Williams (WW) number density or flux of photons from the Pb nucleus as a function of the photon energy k ± = (M Υ /2) exp(±y) with M Υ being the mass of the Υ meson.The cross sections for the underlying photoproduction subprocesses are labelled by σ Pbγ (k − )→PbΥ and σ γ (k + )Pb→Υ Pb .The two terms in Eq. ( 1) correspond to the right-moving and left-moving photon sources, which results in a two-fold ambiguity of the photon energy at a given value of y = 0.
The WW flux is given by a convolution of the impactparameter dependent photon flux N A γ (k, b) calculable in QED [34]

and the nuclear suppression factor
Here, b is the two-dimensional vector between the centres of the two colliding Pb nuclei in the transverse plane, and A A ( b) encodes the Glauber-model probability of having no additional hadronic interaction in the event; for details see [24].
The cross section for the photoproduction process mediating the ultraperipheral Pb + Pb → Pb + Υ + Pb reaction can be expressed in terms of the exclusive photoproduction cross section per bound nucleon N , dσ γ N →Υ N A (W )/dt, and the nuclear form factor F A (t) as where is the t-differential cross section evaluated at t = 0, the variable t being the squared momentum transfer in the process, W is the γ -N c.m.s.energy, and 2 is the minimal momentum transfer squared with m N denoting the nucleon mass.
The nuclear form factor F A (t) is well known from measurements of elastic electron-nucleus scattering and for heavy nuclei it is typically given by the Fourier transform of the two-parameter Woods-Saxon charge distribution ρ(r ) [35], where with |q| = √ −t.We take d = 0.546 fm for the nucleus skin depth and R A /fm = 1.12A 1/3 − 0.86A −1/3 for the nuclear radius.The normalization ρ 0 ≈ 0.17 fm −3 is fixed by requiring that F A (0) = A = 208 for Pb.
The hard scattering amplitude for exclusive Υ photoproduction per nucleon N bound in the nucleus A can be described at NLO in collinear factorization by [6], where In Eq. ( 7), e b = 1/3 and m b = M Υ /2 are the electric charge and the mass of the bottom quark, respectively; α is the fine-structure constant; N c = 3 is the number of colors; γ and * Υ are the polarization vectors of the initial-state photon and the final-state vector meson, respectively; O 1 Υ is the non-relativistic QCD (NRQCD) matrix element for the Υ → b b transition, which is proportional to the radial Υ wavefunction at the origin and which is fixed by the experimental value of the Υ decay width to a dilepton pair, see [36].Note that in this approach, M Υ = 2m b .
The reduced matrix element I (ξ, t = 0) is given by a convolution of the gluon T g (x, ξ, μ R , μ F ) and quark T q (x, ξ, μ R , μ F ) NLO coefficient functions with the gluon F g (x, ξ, t, μ F ) and quark singlet F q,S (x, ξ, t, μ F ) matrix elements involving the corresponding GPDs.Note that the coefficient functions depend on the longitudinal momentum fraction x, the skewness ξ = M 2 Υ /(2W 2 − M 2 Υ ), the renormalization scale μ R , and the factorization scale μ F .In our analysis, we set μ = μ R = μ F and vary μ in the m b /2 ≤ μ ≤ 2m b interval.
In the leading-twist approximation and neglecting the mass of the nucleons, the factors F g and F q,S in the t = 0 limit can be expressed in terms of the helicity-conserving gluon H g (x, ξ, t, μ F ) and quark singlet H q,S (x, ξ, t, μ F ) GPDs as follows [10], with At ξ = t = 0, these GPDs reduce to the usual gluon, quark, and antiquark PDFs of the (bound) nucleons, where x ∈ [0, 1].In Eq. ( 8), each value of the skewness parameter ξ entails an integration over the convolution variable x.In the literature, see [10] for review, the |x| ≥ ξ interval is called the DGLAP region since GPDs there can be interpreted as parton distribution functions evolving in log(μ 2 F ) according to the modified Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations.The |x| < ξ interval is called the ERBL region because GPDs there resemble parton distribution amplitudes, whose μ F evolution is given by the modified Efremov-Radyushkin-Brodsky-Lepage (ERBL) evolution equations.In this work, we employ the Shuvaev transform at NLO to model the ξ dependence of GPDs and to relate the GPDs to PDFs in the DGLAP region, see details in Sect.2.2.To counteract the possible invalidity of the Shuvaev transform in the time-like ERBL region of |x| < ξ, we convolute the GPDs with only the imaginary part of the gluon and quark coefficient functions in Eq. ( 8), which vanish identically for |x| < ξ.We then restore the real part via the high-energy dispersion relation [37] We have checked that this relation accurately reproduces the directly computed real part contribution for W 40 GeV at a percent level in the case that GPDs are approximated by PDFs.At smaller W the deviation increases but, as will be discussed, our main data-driven predictions will nevertheless be valid only for W 100 GeV.
To summarize, the standard pQCD approach to the calculation of the σ γ A→Υ A (W ) cross section is based on Eqs. ( 3) and ( 4), where the hard scattering nuclear amplitude per bound nucleon M γ N →Υ N A (ξ, t = 0) is calculated using the bound nucleon (nucleus) gluon and quark GPDs, see Eqs. (7) and (8).Replacing the bound nucleon by the free proton in these equations, one readily obtains the NLO pQCD predictions for the proton target.The cross section of exclusive Υ photoproduction on the proton reads [compare to Eq. ( 3)] where B Υ (W ) is the energy-dependent slope of the t dependence of the γ + p → Υ + p cross section, which is assumed to be exponential; dσ γ p→Υ p (W )/dt (t = 0) is the differential cross section at t = 0, which is calculated using Eqs.( 4), ( 7) and ( 8) with nuclear GPDs replaced by their free-proton counterparts.
The t dependence of the γ + p → Υ + p cross section has never been measured.Therefore, for the B Υ (W ) slope, we use the following parametrization motivated by Regge phenomenology, where B 0 = 4.63 GeV −2 , α IP = 0.06 GeV −2 , and W 0 = 90 GeV.While the value of B 0 is compatible with fits to the t dependence of elastic J/ψ photoproduction on the proton at HERA [26,38], the value of slope of the Pomeron trajectory α IP is fixed by Model 4 of [39], which fits a wide variety of data on diffraction in proton-proton scattering at the LHC.

GPD modeling
Generalized parton distributions naturally appear in the framework of collinear factorization for hard exclusive processes [7] and combine properties of usual PDFs, distribution amplitudes and elastic form factors [8][9][10].Since GPDs depend on two light-cone momentum fractions x and ξ , the invariant momentum transfer squared t, and the factorization scale μ F , their modeling and extraction from the available experimental data has been notoriously challenging, see, e.g.[40].However, at small values of the skewness ξ , GPDs rather closely resemble usual PDFs in the |x| ≥ ξ DGLAP region and the |x| < ξ ERBL region plays typically only a minor role.These facts significantly simplify the modeling of GPD-originating effects even if the experimental constraints for the three-dimensional structure of GPDs are weak.One of the most widely used models of GPDs at small ξ is based on the so-called Shuvaev transform, which is a method to analytically solve the LO Q 2 evolution equations of GPDs [31][32][33].It is a generalization of solving the usual DGLAP evolution equations using Mellin moments of PDFs.To briefly summarize the method, one first defines effective PDFs, whose Mellin moments are equal to the Gegenbauer (conformal) moments of GPDs.One then inverts these relations and expresses GPDs as certain integrals of the effective PDFs at any given factorization scale μ F .Finally, using the condition of polynomiality of the conformal moments (see details in [32,41]), one argues that the effective PDFs can be approximated by the usual PDFs and obtains the desired connection between GPDs at small ξ and PDFs.In other words, the input GPDs at some low scale μ 0 are assumed to be independent of ξ , and the ξ dependence is then generated radiatively during the scale evolution -this warrants to speak about perturbative skewness.Moreover, since the mixing of the conformal moments under the NLO Q 2 evolution is suppressed by powers of ξ , the Shuvaev transform can also be safely used at NLO in the ξ 1 limit [31].As a phenomenological application of the method, it was shown in NLO and next-to-next-to-leading order (NNLO) analyses [42] that a flexible parametrization of quark and gluon GPDs of the proton in terms of their conformal moments describes well the available HERA data on deeply virtual Compton scattering (DVCS) on the proton.In the case that the condition ξ 1 is not met, the Shuvaev transform should be substituted by explicitly solving the GPD evolution equations [43,44].
One should also point out several general theoretical limitations of the Shuvaev transform.First, approximating effective PDFs by usual PDFs manifestly violates the support area [45], which however is expected to be numerically small at small ξ .Second, to obtain GPDs in the momentumfraction space from the conformal moments, the moments must be analytically continued into the complex Mellin N plane.If there exist singularities in the right half of this plane (Re(N ) > 1), the inversion to obtain x-space GPDs from the conformal moments may prove troublesome or be less accurate [46].This can in principle happen in both DGLAP and ERBL regions.However, in the case of the DGLAP region, it was argued in Ref. [41] that singularities in the right half cannot arise on grounds of Regge theory and the form of the input distribution.To circumvent the potential problem in the ERBL region, we do not use the Shuvaev transform there and instead restore the real part of the amplitude through the dispersion relation, see Eq. ( 12) and its discussion above.
In our work, we employ the Shuvaev transform at NLO as a means to relate the GPDs to PDFs in the DGLAP region.Thus, the quark and gluon GPDs are obtained as integrals of the corresponding quark and gluon PDFs, where the kernel of the transform is As we explained above, Eq. ( 15) is used only to calculate the imaginary part of the hard scattering amplitude (ξ, t = 0).The real part probing the ERBL region is restored via the high-energy dispersion relation (12).In practice, the Shuvaev integrals in Eq. ( 15) involving derivatives of the input PDFs converge rather slowly and have to be precomputed before evaluating Eq. (8).To this end, we have computed the GPDs in a three-dimensional x, ξ/x, μ 2 grid using Eq.(15).The construction of the GPD grid is optimised such that areas in the parameter space that result in a flat interpolation are not overly populated: having more points around ξ/x ∼ 1 mitigates edge effects at the boundary of the DGLAP and ERBL regions [41], while the interpolation in μ 2 is relatively smooth and requires fewer points.
In Fig. 1, we illustrate the effect of finite skewness in GPDs by comparing the gluon and quark-singlet GPDs, F g (x, ξ) and F q,S (x, ξ), obtained through the Shuvaev transform, with their values at ξ = 0, xg(x) and q S (x), as a function of x at the scale μ F = m b .We have used here the CT18ANLO proton PDFs [47] taken from the LHAPDF library [48].In these plots, we have fixed ξ ≈ 10 −3 , which corresponds to the kinematic value of the skewness parameter probed in Υ photoproduction in Pb + Pb UPCs at 5.02 TeV and y = 0.The distributions are plotted in a small interval of the DGLAP region, x ∈ [ξ, 10 −2 ], where the Shuvaev transform is a reliable way to obtain the perturbatively generated skewness of the GPDs.One can see from the figure that the effect of skewness -the deviation between the blue and orange curvesis rather small for most values of x, but grows towards the point x = ξ , especially in the case of quarks.At the same time, to compare with the commonly used skewness factor due to the Shuvaev transform [32,49], we also show the gluon and quark singlet PDFs evaluated at the x + ξ point, (x +ξ)g(x +ξ, μ F ) and q S (x +ξ, μ F ).In this case, the effect of skewness is noticeable (the deviation between the blue and green lines is significant).However, in our NLO pQCD analysis the (x + ξ)g(x + ξ, μ F ) and q S (x + ξ, μ F ) PDFs do not play any special role and we find that the numerical effect of the skewness effects induced by the Shuvaev transform in the calculated cross sections of Υ photoproduction on the proton and a heavy nucleus is small.
The effect of the Shuvaev transform is larger at larger scales μ F , where the effective power growth of the partons becomes steeper and reflects the sensitivity of the Shuvaev transform to the slope of the input PDFs through Eq. (15).The enhancement in the quark singlet GPD is clearly larger than that in the gluon one.However, our analysis shows that the contribution of the quarks is subleading and so the overall effect of incorporating the skewness through the Shuvaev transform is dictated by the gluon GPDs.Our results for the differences between F g (x, ξ, μ F ) and xg(x, μ F ), and F q,S (x, ξ, μ F ) and q S (x, μ F ), are qualitatively similar to those presented in the DGLAP region at LO in Fig. 3 of Ref. [43].
To further support our conclusion that the effects due to GPD modeling turn out to be small in our case, we also examined the gluon and quark singlet GPDs using the double distribution (DD) model of GPDs [9], which has the correct forward limit, satisfies the property of polynomiality (the issue of the D-term does not affect the DGLAP region), and leads to fruitful phenomenology [40,50].Employing the CT18ANLO PDFs in the standard DD profile function [50] at Q 0 = 1.5 GeV, we performed the Q 2 evolution of the GPDs [49,51] to the scale of interest μ F = m b , and confirmed that the predictions of the DD model of GPDs and the Fig. 1 The gluon (left panel) and quark singlet (right panel) GPDs (blue curves) F g (x, ξ) and F q,S (x, ξ) with ξ ≈ 10 −3 obtained through the Shuvaev transformation, compared with PDFs xg(x) and q S (x) (orange dashed curves) at μ F = m b as a function of x.In addition, we also present the distributions (x + ξ)g(x + ξ) and q S (x + ξ) (green dotted curves) Shuvaev transform are indeed rather close.This agrees with a general observation that the Q 2 evolution of GPDs washes out information on their shape at the input scale [44,49], which is also supported by an analysis within the framework of the flexible parametrization of GPDs in terms of their conformal moments [42].

Data-driven approach
As discussed in our previous works in the context of exclusive J/ψ photoproduction in Pb + Pb and O + O UPCs [24,25], the photoproduction scattering amplitude M γ N →Υ N A (ξ, t = 0) introduced above suffers from a large factorization/renormalization scale dependence.While it is milder for the case of Υ photoproduction considered here since the interaction scale is higher than in the J/ψ photoproduction, it is still rather sizeable as we will show later on in Sect. 4. In addition, the NLO results will be shown to somewhat underpredict the HERA and LHC data on the γ + p → Υ + p cross section.An approach to alleviate the strong scale dependence through consideration of additional power corrections ∼ O(μ 2 F /Q 2 0 ) arising in the so-called Q 0 subtraction, where Q 0 is the PDF or GPD parametrization scale, was advocated in [52][53][54][55][56] in the context of p + p and p + Pb collisions.Instead of the Q 0 subtraction, we adopt a data-driven pQCD approach, where the γ + Pb → Υ + Pb cross section is given by the product of the ratio between the Υ photoproduction cross sections on the nucleus and the proton calculated in NLO in pQCD, and the γ + p → Υ + p cross section fitted to the available HERA [26][27][28] and LHC data [29], Using a simple power-like ansatz for σ γ p→Υ p fit (W ) with an additional factor parametrizing the behavior of the cross section near the kinematic threshold [57], one obtains [58] with W 0 = 100 GeV.Note that while the 2018 CMS data [30] have not been included in the fit, they are nevertheless well reproduced, see Fig. 2 ahead.One way to interpret Eq. ( 17) is that we supplement the fitted γ + p → Υ + p cross sections by the theoretical nuclear modification R(W ), which can be anticipated to carry a reduced dependence on the choice of the factorization scale and on the explicit modeling of GPDs.In the first approximation, these effects cancel in R(W ) and it becomes mainly sensitive to the PDFs of protons and nuclei.Alternatively, one can interpret that in Eq. ( 17) one rescales the calculated γ + Pb → Υ + Pb cross sections by a factor that is needed to match the calculated γ + p → Υ + p cross sections with the experimental onesan effective "K factor".In what follows, we will call the cross sections computed through Eq. ( 17) the "data-driven" ones, in contrast to the "standard" pQCD predictions calculated without any reference to experimental data.The approach here is similar in spirit to the leading-order pQCD analysis of the nuclear suppression factor for exclusive J/ψ photoproduction in Pb + Pb collisions introduced and discussed in Refs.[57,59,60].

Results
In this section, we present and discuss our results for the Υ photoproduction process on the proton, γ + p → Υ + p, and the rapidity-differential Υ spectra in Pb + Pb UPCs, Pb + Pb → Pb + Υ + Pb.To estimate the sensitivity of our predictions to higher-order perturbative corrections, we adopt a standard, conservative prescription and vary the factorization and renormalization scales together in the interval of μ F = μ R ∈ {1/2, 1, 2}×m b .As input proton and nuclear PDFs, we use CT18ANLO [47] and EPPS21 [21] PDFs, respectively, from the LHAPDF interface [48].The corresponding GPDs are obtained using the Shuvaev transform as discussed in Sect.2.2.Note that we use the version "A" of the CT18NLO analysis since this was the free proton baseline used in the EPPS21 nPDF analysis.It differs from the default CT18NLO mainly in the strange quark distributions.
In the first instance we make NLO predictions following the standard pQCD approach, and then subsequently compare and contrast features of these predictions with those obtained from the data-driven method explained in Sect.3, as well as with our earlier analyses [24,25] of J/ψ photoproduction in Pb + Pb UPCs.

Standard pQCD results for γ + p → Υ + p cross section
Figure 2 presents the σ γ p→Υ p (W ) cross section of exclusive Υ photoproduction on the proton, γ + p → Υ + p, as a function of the invariant photon-proton c.m.s.energy W .The dashed, dot-dashed and dotted curves correspond to the NLO pQCD predictions of Eq. ( 13), which as input use either the proton GPDs obtained via the Shuvaev transform (the curves labeled "GPD") or the usual proton PDFs, i.e., the ξ = 0 forward limit of the GPDs (the curves labeled "PDF").Each pair of predictions is evaluated with three scale settings μ = μ F = μ R ∈ {1/2, 1, 2} × m b .The shaded band represents the propagated uncertainty of the proton PDFs used for the GPD-based predictions at μ = m b .These results are compared with the available HERA [26][27][28] and the LHC data [29,30] on this process.Note that it is argued in [55] that the extracted values of σ γ p→Υ p (W ) at the largest W from the LHCb rapidity-differential measurements [29] should be shifted upwards because the collaboration used a less accurate approximation for the photon flux in their analysis.Finally, the black solid line labeled "Fit" is the parametrization of Eq. (18).One can see from the figure that while our NLO pQCD predictions reproduce the trends of the W dependence of the data, they underestimate the normalization of the cross section, especially at larger values of μ.A reliable description of the normalization would thus require a better theoretical understanding of the perturbative structure of the process including, e.g. the relevance of unknown next-to-NLO corrections, significance of the double logarithmic α s log(μ 2 F /m 2 b ) log(1/ξ ) terms present already in the NLO hard coefficient functions T g and T q [6,61], 1 and the size of the relativistic corrections to the quarkonium wave function [62].An account of these effects is beyond the scope of this paper and we will work around these issues through the data-driven predictions.
The systematics of the NLO pQCD predictions in Fig. 2 can be summarized as follows.First, as discussed in Sect.2.2, the effect of skewness is rather mild, i.e., the difference between the GPD-based and PDF-based predictions is small, especially at smaller values of μ.Second, while the GPD-based predictions correspond to higher values of σ γ p→Υ p (W ) than the corresponding PDF-based ones at μ = m b and μ = 2m b , this hierarchy of predictions is reversed at μ = m b /2.A detailed examination indicates that this originates from a delicate interplay among the LO gluon and NLO gluon and quark contributions in M γ N →Υ N A (ξ, t = 0) whose relative signs vary depending on the scale choices.This is further complicated by the fact that the magnitude of the skewness effect generated by the Shuvaev transform ( 15) depends on both W (through its dependence on ξ ) and μ F controlling the slope of the x dependence of the gluon and quark PDFs.Third, as a result of scale-dependent sign differences of quark/gluon contributions, the relative ordering of predictions from low to high μ depends on W .  3), ( 4) and (7).As input, we use the nuclear GPDs constructed using the Shuvaev transform and the EPPS21 nPDFs (central plus error sets).The three curves correspond to the three different choices of the factorization/renormalization scales μ = {m b /2, m b , 2m b }.The shaded band gives the propagated uncertainty of the EPPS21 nPDFs in the μ = m b case.As a useful reference, the upper x-axis shows the values of W + corresponding to each y, that is, W + = (M Υ √ s N N e y ) 1/2 .Two features of the results in Fig. 3 deserve to be mentioned.First, one can see from the figure that apart from the very tails of the rapidity distribution, |y| > 3, the central prediction with μ = m b does not lie between the other scale choice predictions with μ = m b /2 and μ = 2m b .This feature can be readily observed also in the results for the proton cross section in Fig. 2. Indeed, taking, for instance, The HERA [26][27][28] and LHC data [29,30] data on this process are shown as well, together with a fit [Eq.( 18)] to the HERA data (the black solid line labeled "Fit") y = 0 corresponding to W ≈ 200 GeV, one can see that the predictions for σ γ p→Υ p (W ) with μ = m b lie below the corresponding predictions at μ F = m b /2 and μ F = 2m b .Second, the scale uncertainty is rather large and the prediction with μ = m b /2 lies outside the nPDF uncertainty band.We will show in Sect.4.3 that both of these features can be tamed through our data-driven approach.
In Figs. 4, 5 and 6, we show various decompositions of dσ Pb+Pb→Pb+Υ +Pb /dy at μ = m b as a function of y. Figure 4 presents the breakdown of the full cross section into the quark, gluon and interference contributions.It is clear that over the entire considered rapidity region, the gluon contribution dominates the quark contribution, in dissimilarity to the analogous breakdown for the J/ψ rapidity-differential cross section in Pb + Pb UPCs in NLO pQCD shown in our previous studies [24,25], where the quark contribution was shown to be the dominant one around mid rapidity.One should note that even if the quark contribution is small, it is not zero or structureless and it leads to a visible contribution in the interference terms.One can speculate that the interaction scale in the Υ photoproduction is already sufficiently large so that NNLO corrections will not change the mutual hierarchy of quark/gluon contributions.The situation could be very different in the case of J/ψ photoproduction where the quark dominance is a consequence of a coincidental cancellation between the LO and NLO gluon contributions.
In Fig. 5, we show the W + and W − decomposition i.e. separately plot the two contributions in Eq. ( 1).The situation is very similar to that in our J/ψ analysis, see [24,25] for more details.For instance, the W − contribution dominates at positive forward rapidities (large W + ) because there (kd N Pb γ /dk) k=k − (kd N Pb γ /dk) k=k + .The situation is reversed in the region of backward rapidities corresponding to small W + .The presence of two terms in Eq. ( 1) complicates the extraction of the small-x contribution from UPC cross sections at y = 0.However, it is possible to separate the W + and W − contributions by studying UPCs accompanied by forward neutron emission due to electromagnetic excitation of one or both colliding nuclei [63].Such an analysis in the case of coherent J/ψ photoproduction in Pb + Pb UPCs at 5.02 TeV was recently performed by the CMS collaboration [64], which allowed one to deepen the small-x reach down to x ∼ 10 −4 .
Finally, Fig. 6 presents the decomposition of dσ Pb+Pb→Pb+Υ +Pb /dy into the contributions of the real and imaginary parts of M γ N →Υ N A (ξ, t = 0).The imaginary part clearly dominates over the entire range of rapidity.Again, the situation was much more involved in the case of J/ψ, where the interplay of the two was highly non-trivial [24,25].

Data-driven pQCD predictions for the
Pb + Pb → Pb + Υ + Pb UPC cross section The data-driven pQCD prediction for the UPC cross section dσ Pb+Pb→Pb+Υ +Pb /dy is given by Eq. ( 17), where only the ratio of the nucleus and proton cross sections R(W ) of Eq. ( 19) is calculated using our NLO pQCD framework, while the absolute normalization is given by σ γ p→Υ p fit (W ) obtained from a fit to the proton data, see Eq. ( 18).The results for the differential cross section as a function of the Υ rapidity y are shown in Fig. 7.The numerator and the denominator of the ratio R(W ) are calculated using the EPPS21-based and the CT18ANLO-based GPDs, respectively; these curves are labeled "nGPD".For comparison, we also show the results of the calculation, where we neglect the effect of skewness and use the forward ξ → 0 limit for nuclear and proton GPDs; these curves are labeled "nPDF".The blue dot-dashed curve represents our central prediction at μ = m b with the blue shaded band quantifying the propagation of the EPPS21 nPDF and the CT18ANLO proton PDF uncertainties; their counterparts in the case, where GPDs are taken in the forward limit, are given by the red solid curve and the corresponding red shaded band.The dotted and dashed curves correspond to the ratio R(W ) evaluated at μ = m b /2 and (W ) fit to the γ + p → Υ + p photoproduction data is an extrapolation: the HERA data are available only for W ≥ 100 GeV (see Fig. 2), but for |y| ≥ 2 there is a large contribution from W < 100 GeV, see Fig. 5.
It is important to contrast our results in Fig. 7 with the standard NLO pQCD predictions shown in Fig. 3. First, while the shapes of the y distribution are very similar, the normalization of the data-driven results is approximately a factor of 2-2.5 higher.This is a straightforward consequence of the rescaling of the cross section of exclusive Υ photoproduction on the proton to fit the available data.Second, the dependence on the factorization/renormalization scale μ is now more regular in the central rapidities: the central prediction with μ = m b lies below the μ = m b /2 result and above the μ = 2m b one.Most importantly, the scale dependence has reduced significantly.Third, the effects of GPD modeling are seen to largely cancel in the ratio R(W ).As a result, the data-driven pQCD predictions for dσ Pb+Pb→Pb+Υ +Pb /dy are here mainly sensitive to the input PDFs.Note that for lower W , where the real part restoration via Eq.( 12) is less accurate, the behavior of our results is less regular, but this lies in the tails of the y distributions where our predictions anyhow lean on an extrapolation of σ γ p→Υ p fit (W ) into non-measured values of W .
To quantify the magnitude of nuclear effects probed in exclusive Υ photoproduction in Pb + Pb UPCs at the LHC, it is convenient to consider separately the ratio R(W ) in Eq. (19).Indeed, at a given value of the Υ rapidity y = 0, the Pb + Pb UPC cross section contains two terms leading to a two-fold ambiguity in the photon-nucleon c.m.s.energy W ± .As a consequence, this mixes the low-x and medium-x contributions to the UPC cross section and makes it challenging to extract the information on small-x physics, which is often thought to be at the heart of the process under consideration.This issue is absent in the case of R(W ) although it cannot be experimentally measured in a model-independent way.In the upper panel of Fig. 8, we show the ratio R(W ) as a function of W + .On the x-axis at the top, we also give the corresponding values of the skewness The curve corresponds to the central prediction at μ = m b shown in Fig. 7, where the numerator and the denominator of R(W ) are calculated using the EPP21-based nuclear GPDs and the CT18ANLO-based free proton GPDs, respectively.The shaded band is the result of the propagation of the EPPS21 nPDF and the CT18ANLO proton PDF uncertainties.We see that the rescaling factor R(W ) depends strongly on W and its value can be as large as several hundreds.The absolute value can, however, be mostly explained through the proton and nuclear form factors.To see this and to provide a closer comparison with nuclear modifications of nPDFs, one can eliminate the effects of the nuclear and the proton form factors in the R(W ) ratio by rescaling it by the factor of R (W ), where B Υ (W ) is the slope of the t dependence of the γ + p → Υ + p differential cross section in Eq. ( 14) and F A (t) is the nuclear form factor in Eq. ( 5).Note that R (W ) depends on W + through |t min | = m 2 N (M Υ /W + ) 4 and B Υ (W + ).In the lower panel of Fig. 8, we present the scaled R(W ) ratio, i.e., the product R(W )× R (W ), as a function of W + by the red solid curve.The propagated nuclear and free proton PDF uncertainties are given by the red shaded band.One can see from the figure that as a function of ξ + , R(W ) × R (W ) exhibits significant suppression for small ξ + < 0.05 and a ∼ 10% enhancement at ξ + ∼ 0.1.This behaviour reflects the characteristic nuclear modifications of nPDFs associated with nuclear shadowing at small x and as a function of ξ = ξ + , where g A (q S A ) and g p (q S p ) are the gluon (quark-singlet) distributions per nucleon in the nucleus and the free proton, respectively.The corresponding shaded bands represent the EPPS21 nPDF uncertainties of these ratios.One can see that the shape and normalization of both R 2 g (ξ ) and R 2 q (ξ ) is similar to those of R(W ) × R (W ).Moreover, because of the dominance of the gluon-initiated contribution over the quark one, see Fig. 4, and the flat shape of the gluon nuclear modifications at small x, the values of R(W )× R (W ) and R 2 g (ξ ) become very close for ξ + ≤ 10 −3 (W + > 200 GeV).

Feasibility of the measurement of Υ photoproduction
in Pb + Pb UPCs at the LHC Having now obtained an educated estimate for the Υ cross section in Pb + Pb collisions, we will here check to what extent an experimental measurement of the process would be feasible.To this end, we lean on the exclusive Υ p + Pb measurement by the CMS collaboration [30] at √ s N N = 5.02 TeV.This measurement with an integrated luminosity of L( p + Pb) = 32.6 nb −1 reported ∼ 80 identified Υ (1S) particles and yielded a total cross section σ ( p + Pb) = 94.8 nb in the rapidity interval |y| < 2.2.If we desire a Pb + Pb measurement that is as precise as the p + Pb measurement (the same number of events), and assume the same efficiency, the condition is From Fig. 7, we find a total cross section σ (Pb+Pb) ∼ 52 µb in the same rapidity interval −2.2 < y < 2.2.It then follows that the required integrated luminosity should be to observe ∼ 80 events.Given that the recorded luminosity at the 2018 Pb + Pb run for CMS is as high as 1.7 nb −1 [65], our counting would thus promise ∼ 80 × (1.7/0.06)≈ 2300 events.Moreover, in Run III CMS aims for an integrated luminosity of 13 nb −1 [19], so the measurement of exclusive Υ photoproduction in Pb + Pb collisions looks more than feasible to be performed at the CMS experiment.

Conclusions and outlook
We presented the first study of the rapidity-differential cross section of exclusive Υ photoproduction in ultraperipheral lead-lead collisions at the LHC using collinear factorization at NLO pQCD.In addition, we extended our previous framework in [24,25] by now including explicit GPD modeling through the Shuvaev integral transform.In our standard NLO pQCD approach, we showed that the GPD effects are small, and unlike in the J/ψ case, the imaginary part and gluon contributions dominate the amplitude.The scale uncertainties are significantly reduced from the J/ψ case, but they are still alarmingly large.In the γ + p case, the NLO calculation was shown to underpredict the HERA data, which calls for further improvements such as NNLO pQCD and NRQCD corrections.
Using the γ + p → Υ + p cross section from HERA and the LHC as a baseline, we proposed a data-driven pQCD approach to make more constrained predictions for dσ (Pb + Pb → Pb + Υ + Pb)/dy and showed that the resulting factorization/renormalization dependence becomes smaller than that in the standard pQCD result for this process.In addition, effects due to the explicit GPD modeling largely cancel and most of the remaining uncertainty is due to PDFs of free and bound nucleons.This serves as a first step towards being able to include heavy quarkonia UPC data in the global analyses of nPDFs to provide constraints on partons inside nuclei at moderate to low x.We also estimated that the production cross sections are high enough for this process to be measured in Pb + Pb collisions at the LHC.While the theoretical situation nevertheless seems a little better for the Υ production, the experimental statistics obtainable may be sparser than that for J/ψ.Future works can therefore include applying our data-driven approach to exclusive J/ψ photoproduction in nucleus-nucleus collisions, where also the statistical quality of the baseline γ + p data is greater than for Υ production.In the J/ψ case, the GPD modeling given by the Shuvaev transform is surmised to have an even smaller effect in comparison to the Υ photoproduction considered here, but to what extent the scale dependence can be tamed, calls for a detailed analysis.
Data Availability Statement This manuscript has no associated data or the data will not be deposited.[Authors' comment: The analysis output data are available from the authors upon request.]Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.Funded by SCOAP 3 .SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.

Fig. 2
Fig.2The γ + p → Υ + p cross section as a function of W .The NLO pQCD GPD-based (red curves) and PDF-based (blue curves) predictions evaluated at μ = {m b /2, m b , 2m b } are presented by the dashed, dot-dashed and dotted lines; the shaded band is the propagated CT18ANLO PDF uncertainty for the GPD-based result at μ = m b .The HERA[26][27][28] and LHC data[29,30] data on this process are shown as well, together with a fit [Eq.(18)] to the HERA data (the black solid line labeled "Fit")

Fig. 3 Fig. 4 Fig. 5 Fig. 6
Fig. 3 Standard NLO pQCD prediction for the rapidity-differential cross section for exclusive coherent Υ photoproduction in Pb + Pb UPCs as a function of the Υ rapidity y at √ s N N = 5.02 TeV.We use nuclear GPDs constructed from the EPPS21 nPDFs via the Shuvaev transform.The dashed-dotted curve represents the prediction with

Fig. 7
Fig. 7 Data-driven NLO pQCD prediction for the rapidity-differential cross section for exclusive coherent Υ photoproduction in Pb + Pb UPCs as a function of the Υ rapidity y at √ s N N = 5.02 TeV.We use the nuclear and proton GPDs constructed from the EPPS21 nPDFs and CT18ANLO proton PDFs, respectively, obtained via the Shuvaev transform (the curves labeled "nGPD").For comparison, we also show the results based on the ξ = 0 limit of the used GPDs (the curves labeled

Fig. 8
Fig. 8 Upper panel: The ratio R(W ) = σ γ Pb→Υ Pb (W )/σ γ p→Υ p (W ) pQCD as a function of the c.m.s.energy W + evaluated using the EPPS21 nuclear and CT18ANLO free proton PDFs at μ = m b .The shaded band corresponds to the EPPS21 and CT18ANLO PDFs uncertainties.The upper x-axis indicates the corresponding values of the skewness ξ + .Lower panel: the rescaled ratio R (W ) × R(W ) as a function of W + .For comparison, the EPPS21 gluon and quark-singlet nuclear modifications squared along with their uncertainties are overlaid.The shaded bands show the PDF-originated uncertainties