Next-to-leading order perturbative QCD predictions for exclusive J/ψ photoproduction in oxygen-oxygen and lead-lead collisions at the LHC

,


I. INTRODUCTION
Traditionally parton distribution functions (PDFs) and their nuclear counterparts, nuclear PDFs (nPDFs), have been determined from inclusive processes such as lepton-hadron deep inelastic scattering (DIS) and the production of leptons pairs (Drell-Yan process), light and heavy mesons, dijets, and gauge bosons in hadron-hadron scattering, see Refs.[1][2][3][4] for recent reviews.The determination of proton and nuclear PDFs has become an active branch of phenomenological applications of quantum chromodynamics (QCD), for recent examples of global fits of PDFs and nPDFs, see [5][6][7][8][9][10][11][12][13].However, despite the dramatic progress in the methodology of PDF extraction from the available data, including an account of higherorder (up to next-to-next-to-leading order, NNLO) perturbative QCD corrections, effects of heavy (charm and bottom) quark masses and small-x resummation and the reliance on sophisticated statistical and computational methods (Bayesian and Hessian error estimates and neural networks), the resulting PDFs and nPDFs still suffer from significant uncertainties.
As a consequence, there is a continuing interest to explore novel kinematics, processes, and observables, which would allow one to obtain additional constraints on the PDFs.In particular, it has been discussed that the exclusive photoproduction of J/ψ mesons on the proton and nuclear targets in the so-called ultraperipheral collisions (UPCs) allows one to probe the gluon density of the target at small momentum fractions x ∼ 10 −5 − 10 −3 and resolution scales µ 2 ∼ 3 GeV 2 [14][15][16][17][18][19] (photoproduction of other quarkonium states, ψ ′ and Υ, has also been considered).This is based on the early observation that in the leading logarithmic approximation, i.e., to the leading order (LO) of perturbative QCD (pQCD), the cross section of this process is directly proportional to the gluon density squared [20].However, it was later found that the next-to-leading order (NLO) QCD corrections involving both gluon and quark distributions are very large [15,21], which questions the common interpretation in terms of the gluon density.While several methods to stabilize the NLO results have been proposed [16,22,23], further theoretical and phenomenological studies are still required.
We recently performed a detailed study of the cross section of exclusive photoproduction of J/ψ mesons in Pb-Pb UPCs at the Large Hadron Collider (LHC) as a function of the J/ψ rapidity y in the framework of collinear factorization and NLO pQCD, and confirmed the dramatic role of the NLO effects [24].In particular, we found that at central rapidities the cross section is dominated by the quark contribution since the gluon one largely cancels in the sum of the LO and the NLO terms.Additionally, even though the scale dependence of our results turned out to be -as expected -rather sizable, we determined an "optimal scale" allowing for a simultaneous reasonable description of the available Run 1 and Run 2 LHC data on this process.In addition, we observed that the amplitude for this process is predominantly imaginary in a broad range of rapidities with a small window at forward and backward rapidities, where the real part gives the dominant contribution.
The purpose of this work is to extend the analysis of Ref. [24] by (i) updating our previous results for Pb-Pb collisions with three different state-of-the-art nPDF sets, namely, EPPS21 [11], nNNPDF3.0[12], and nCTEQ15WZSIH [8], (ii) making detailed predictions for the O + O → O + J/ψ + O rapidity-differential cross section for the planned oxygen run at the LHC [25,26], and (iii) presenting predictions for the ratio of the J/ψ rapidity distributions in Pb-Pb and O-O UPCs.This allows us not only to better control the theoretical uncertainties associated with the nPDFs, but also to tame (reduce) the scale dependence of our NLO results by considering the ratio of J/ψ production with different nuclear collision systems.For the recent predictions of J/ψ photoproduction in O-O UPCs at the LHC in the color dipole framework, see [27].
The rest of the paper is organized as follows.In Sec.II, we recapitulate the framework of NLO pQCD coherent exclusive photoproduction of J/ψ in nucleus-nucleus UPCs, pointing out specific extensions to the oxygen beams.Section III contains our results, which include updated predictions for dσ(Pb+Pb → Pb+J/ψ+Pb)/dy with the most recent sets of nPDFs and their comparison to all available LHC data on this process, detailed predictions for dσ(O + O → O + J/ψ + O)/dy for the oxygen run with an analysis of the scale dependence and the decomposition into the imaginary and real parts as well as into the gluon and quark contributions, and, finally, predictions for the ratios of the J/ψ rapidity distributions in O-O and Pb-Pb UPCs with an exhaustive analysis of the scale and energy dependence.We discuss and summarize our findings in Sec.IV.

II. COHERENT J/ψ PHOTOPRODUCTION IN NUCLEUS-NUCLEUS UPCS IN NLO PQCD
In the equivalent photon approximation the cross section of coherent J/ψ photoproduction in UPCs of nuclei (ions) A 1 and A 2 , as a function of the J/ψ rapidity y, reads [28] where kdN A γ (k)/dk is the flux of equivalent quasi-real photons emitted by ions A 1 and A 2 , k is the photon energy and σ γA→V A (W ) is the cross section of coherent (without nuclear breakup) J/ψ photoproduction on a nuclear target with W being the collision energy of the photon-nucleon system.The two terms in Eq. ( 1) represent two possibilities to arrive at the same final state corresponding either to ion A 1 emitting a photon interacting then with ion A 2 or ion A 2 being a source of photons interacting with ion A 1 .We define the positive rapidity y in the direction of the ion A 1 , from which one obtains that the relation between the photon energies k ± and the rapidity y is k ± = (M J/ψ /2)e ±y , where k + and k − refer to ions A 1 (positive longitudinal momentum) and A 2 (negative longitudinal momentum), respectively, and M J/ψ is the mass of J/ψ.The corresponding photonnucleon system energies are W + = 2M J/ψ e y E 2 and W − = 2M J/ψ e −y E 1 , where E 2 and E 1 are the per nucleon energies of beams A 2 and A 1 , respectively.For symmetric UPCs, we have is the nucleon-nucleon center-of-mass system (c.m.s.) energy.The interference between the amplitudes, where the photons are emitted by different nuclei, is important only at very small values of the momentum transfer t (very small values of the J/ψ transverse momentum) [29] and hence can be safely neglected in the case of the UPC cross section integrated over t which we consider.The flux of equivalent photons emitted by a relativistic ion in UPCs is given by a convolution of the impact parameter dependent photon flux N A γ (k, ⃗ b) and the nuclear suppression factor Γ AA ( ⃗ b), where ⃗ b is a two-dimensional impact parameter vector denoting the distance between the centers of colliding nuclei in the transverse plane.Furthermore, the impact parameter dependent photon flux N A γ (k, ⃗ b) of a relativistic nucleus A with Z protons can be readily calculated in QED [30], where α e.m. is the fine-structure constant, γ L is the nucleus Lorentz factor, J 1 is the cylindrical modified Bessel function of the first kind and FA (t) is the nucleus form factor normalized to one, i.e., FA (t) = F A (t)/A.The nuclear form factor F A (t), accompanied with the normalization condition F A (0) = A, is in turn given by the standard Fourier transform of the nuclear density ρ A (r), where t = −|q| 2 .The nuclear density is well known from measurements of elastic electron-nucleus scattering and is usually parameterized in the form of two-parameter Fermi model (also called Woods-Saxon model) and three-parameter Fermi model (3pF) [31].The former is typical for heavy nuclei, for lead see [24], and latter is usually employed for medium-heavy nuclei.In particular, in this work we use the 3pF parametrization for oxygen with the free parameters determined from nuclear chargedensity measurements [31], d = 0.513 fm and w = −0.051 .
For lead we use d = 0.546 fm and w = 0.The effective nuclear radii are here taken from the following empirical parametrization (see e.g., [32])1 Further, in Eq. ( 2), the nuclear suppression factor Γ AA ( ⃗ b) represents the probability of having no hadronic interactions at impact parameter ⃗ b, which can be evaluated using the Glauber model for nucleus-nucleus scattering where σ N N (s N N ) is the total nucleon-nucleon cross section [33] and For the cross section of coherent J/ψ photoproduction on nuclei A we use the form where the t-dependence governed by the nuclear form factor squared |F A (t)| 2 , is factorized from the cross section of J/ψ production on bound nucleons N of the nuclear target, i.e., dσ γN →J/ψN A /dt (this is indicated by the subscript), where In the case of the tintegrated cross section, the factorized form of Eq. ( 9) approximates with a several-percent precision a more accurate expression that takes into account the correlation between t and x, i.e., the correlation between the momentum transfer and the magnitude of nuclear effects (nuclear shadowing) affecting dσ γN →J/ψN A /dt, see the discussion in [34,35].
The QCD dynamics of the process is contained in the dσ γN →J/ψN A /dt(t = 0) cross section.In the framework of collinear factorization for exclusive processes in NLO perturbative QCD and using the non-relativistic (static) approximation for the charmonium wave function, the cross section reads (see [21,24] for details and references) where the reduced scattering amplitude is given by a convolution of the gluon and quark hard scattering coefficient functions T g (x, ξ, µ R , µ F ) and T q (x, ξ, µ R , µ F ) with the corresponding gluon and quark generalized parton distribution functions (GPDs) F g (x, ξ, t, µ F ) and F q,S (x, ξ, t, µ F ) of the bound nucleons, In Eq. ( 10), e c = 2/3 and m c = M J/ψ /2 is the charm quark mass in the non-relativistic limit, ⟨O 1 ⟩ V is the non-relativistic QCD matrix element associated with the J/ψ → l + l − leptonic decay, ξ = ζ/(2 − ζ) is the socalled skewness parameter with ζ = M 2 J/ψ /W 2 being an analog of Bjorken x in inclusive DIS.Note that the quark contribution in Eq. ( 11) contains a singlet combination of quark GPDs of four active flavors, F q,S (x, ξ, t = 0, µ F ) = q=u,d,s,c F q (x, ξ, t = 0, µ F ).In our analysis we take the factorization and renormalization scales to be equal, i.e. µ = µ F = µ R , which sets the term ∼ β 0 ln(µ 2 R /µ 2 F ) in the NLO gluon contribution to zero, see Eq. (3.72) in [21].We quantify the dependence of our results on the scale choice by varying the scale in the m c ≤ µ ≤ 2m c interval.
In general, GPDs are complicated nonperturbative distributions depending on two light-cone momentum fractions x and ξ and the momentum transfer t as well as the factorization scale µ F .However, in the high-energy limit the skewness parameter is very small (ξ ≪ 1) and its effect on GPDs is expected to be rather moderate.In particular, in the calculation of the dσ γN →J/ψN A /dt(t = 0) cross section, theoretical uncertainties associated with detailed modeling of GPDs are expected to be much smaller than the scale and nPDF uncertainties.Therefore, as a first approximation, we neglect the skewness effect and take GPDs in the forward limit, where they can be identified with the usual PDFs, where now 0 ≤ x ≤ 1.The distributions g A (x, µ F ), q A (x, µ F ) and qA (x, µ F ) are the usual gluon, quark, and antiquark nPDFs per nucleon.In our analysis, we use the recent nPDF sets EPPS21 [11], nNNPDF3.0[12] and nCTEQ15WZSIH [36].

III. RESULTS
A plan is moving forward that an oxygen-oxygen (O-O) run would be performed at the LHC in Run 3 [25,26,37].In addition to shedding light on the soft QCD dynamics and studying hard scattering with small nuclear systems, it might help to address open questions relating to forward scattering physics.From the point of view of UPC studies, as we pointed out in the Introduction, an analysis of coherent J/ψ photoproduction in O-O UPCs and a comparison to the case of Pb-Pb UPCs can be used to constrain theoretical uncertainties of our NLO pQCD predictions for this process.At the time of writing of this paper, it is not yet completely clear at which nucleonnucleon c.m.s.energy √ s N N the O-O run will be completed.Therefore, we will consider four scenarios with √ s N N = 2.76, 5.02, 6.37 [37], and 7 TeV [25], which will help us to better understand the energy dependence of our results.

A. Photon fluxes and nuclear form factors
The results for the kdN A γ (k)/dk photon flux obtained through Eq. ( 2), where k = (M J/ψ /2)e y , i.e. positive rapidity y corresponds to large photon energies k, for oxygen and lead beams for four different values of the invariant collision energy √ s N N are presented in Fig. 1.
In the figure, the blue solid curves correspond to the oxygen case and the orange dashed curves to the lead case.
In order to conveniently compare the two cases, we normalized the fluxes by the factor of 1/Z 2 with Z O = 8 for oxygen and Z Pb = 82 for lead.One can see from the figure that at negative rapidities (small photon momenta) the photon flux of the lead beam is much larger than that for the oxygen beam, At the same time, since the effective nuclear radius of lead is almost 3 times as large as that of oxygen, the spectrum of equivalent photons of lead falls off more rapidly when y is increased (corresponding to an increase in k) than that of oxygen.Eventually, for large values of rapidity y ≥ 4.4 corresponding to k ≥ 120 GeV, the photon flux for oxygen becomes bigger than that of lead.
We have numerically checked that setting w = 0 in Eq. ( 5), i.e., assuming the two-parameter Fermi (2pF) model for oxygen with the same d and R A parameters, leads to a relative difference of under four percent in the photon flux for the photon energies up to k ≈ 50 GeV corresponding to the J/ψ rapidities |y| ≤ 3.5.In addition, we have checked that the photon flux is not sensitive to the exact value of σ N N used in Γ AA ( ⃗ b), see Eq. (8).For example, at √ s N N = 6.37 and 7 TeV, calculations for the photon flux with σ N N = 95 mb and σ N N = 100 mb differ by less than 1% all the way up to |y| = 4. Thus, we conclude that the major difference between the scaled photon fluxes of the oxygen and lead ions originates from the different effective radii R A of these nuclei.
In the left panel of Fig. 2 we show the results for the oxygen and lead form factors scaled by the corresponding mass number A as they are given by Eq. ( 4).The values of the scaled form factors approach one due to our normalization constraint F A (0) = A but, as we move to higher values of t, we see that the scaled oxygen form factor is the dominant one except for the oscillations at very high values of t.But again, since A Pb ≫ A O , the absolute magnitude of the form factor of lead is the bigger one.Then, in the photoproduction cross section σ γA→V A , we have an integral over the square of the absolute value of the form factor.The right panel of Fig. 2 shows the values of this integral scaled by the square of the mass number A for both the oxygen (solid blue) and the lead (dashed orange) cases.Similarly to the photon flux, this ratio gets intertwined in the ratios of the cross sections, but at central rapidities y = 0 corresponding to t min ≈ 10 −6 , we should expect to see a factor of 4.6 from the ratio of the integrals.

B. Rapidity dependent cross sections in Pb-Pb
UPCs and comparison with the LHC data Based on the NLO pQCD theoretical framework outlined in Sec.II, we present below our predictions for the dσ(Pb + Pb → Pb + J/ψ + Pb)/dy cross section of coherent J/ψ photoproduction in Pb-Pb UPCs as a function of the J/ψ rapidity y at √ s N N = 2.76 TeV (Run 1) and √ s N N = 5.02 TeV (Run 2) at the LHC and compare them with all available LHC data on this process.We performed our calculations using the most recent EPPS21 [11], nNNPDF3.0[12], and nCTEQ15WZSIH [36] sets of nPDFs, which updates our predictions in [24].Figure 3 demonstrates the variation of our predictions due to the choice of the scale µ, which is allowed to vary in the m c ≤ µ ≤ 2m c interval (m c = M J/ψ /2 = 1.55 GeV in the nonrelativistic limit that we use): the upper dashed curves correspond to µ = 3.1 GeV, while the lower dotted curves are for µ = 1.55 GeV.The solid curve in each panel corresponds to an "optimal scale", which is chosen to simultaneously describe the central rapidity data available from Run 1 (left panels) and Run 2 (right panels) at the LHC.The Run 1 data at √ s N N = 2.76 TeV include the ALICE data at the central rapidity |y| < 0.9 [38] (labeled "ALICE Cent") and at the forward rapidity 2.6 < |y| < 3.6 [39] labeled "ALICE Forw") as well as the CMS data in the rapidity interval 1.8 < |y| < 2.3 [40] (labeled "CMS Forw").The Run 2 data taken at √ s N N = 5.02 TeV are the ALICE data at midrapidity |y| < 0.8 [41] (labeled "ALICE Cent"), the ALICE data at forward rapidities 2.5 < |y| < 4 [42] (labeled "ALICE Forw"), the LHCb data at forward rapidities 2 < |y| < 4.5 [43] (labeled "LHCb 2015") and their recent update [44] (labeled "LHCb 2018").
A comparison of the results presented in Fig. 3 with our results in [24] shows that the difference between our calculations using EPPS21 and EPPS16 is very small with a very similar value of the optimal scale and the same shape of the y dependence as well as the matching magnitude of the scale dependence and the quality of the data description.To be exact, at central rapidity y = 0, for Run 1 there is a factor of about 22 between the highest scale and the lowest scale results and for Run 2 energy this factor is about 55.
The improvement, when moving from nNNPDF2.0 [45] (Fig. 10 of [24]) to the newer nNNPDF3.0set, is rather dramatic.We find that the shape of the dσ(Pb + Pb → Pb + J/ψ + Pb)/dy cross section at the optimal scale µ = 2.22 GeV is qualitatively similar to that obtained with EPPS16 or EPPS21.Simultaneously, however, the correspondence with the data is slightly worse: while the data at y ≈ 0 is reproduced by construction, the solid curve somewhat underestimates the data at |y| ̸ = 0. Note that the good agreement with the data at y ≈ 0 is im-portant for the comparison of the Pb-Pb and O-O UPC data; see the discussion in Sec.III D. The ratio between the highest scale and the lowest scale at central rapidity is about 17 for Run 1 and about 26 for Run 2.
In contrast to EPPS21 and nNNPDF3.0,we find that the newest nCTEQ15WZSIH nPDF set actually does better on all accounts.The scale dependence at central rapidity is only about a factor of 10 for Run 1 and about a factor of 12 for Run 2. The curve corresponding to the optimal scale of µ = 2.02 GeV goes through the central rapidity data points in addition to the forward/backward data both at Run 1 and Run 2 energies.Moreover, the nCTEQ15WZSIH nPDF set favors the ALICE forward data [42] and the newer LHCb 2018 data [44] over the 2015 LHCb data [43].We have checked that the better agreement of our calculations using the nCTEQ15WZSIH nPDF set with the UPC data is due to the very strongly enhanced strange quark distributions, see Fig. 4 in Ref. [46].Thus, this process may give an interesting opportunity to obtain new constraints on the elusive strange quark distribution in the proton and nuclei.
For all three sets, when considering the full range of scales µ ∈ [m c , M J/ψ ], the scale uncertainty decreases slightly -as was with the earlier EPPS16 set -as we move further away from the central rapidity towards backward and forward rapidities.This is partly because at very large values of rapidity, i.e. |y| > 3, the photoproduction amplitude receives a large contribution from the Wcomponent corresponding to small values of k, which in turn means that we are probing the underlying GPDs at high values of x, where the scale dependence is constrained rather well.In any case, it is interesting to notice that this rapidity dependence seems to be a common property for both the old and the new nPDF sets (see Fig. 4 in [24]).
To estimate the PDF uncertainty of our predictions due to the EPPS21 and nCTEQ15WZSIH nPDFs, we used the following asymmetric form for the uncertainty δO ± [9] where O(S 0 ) denotes the predictions with the central set for the observable O and O(S ± i ) correspond to the values calculated with the plus and minus PDF error sets.In the case of nNNPDF3.0,we used the 90% confidence level (CL) interval prescription [12].All PDF uncertainty calculations are performed at the corresponding values of the optimal scale µ.
Figure 4 illustrates the uncertainty of our predictions for the dσ(Pb + Pb → Pb + J/ψ + Pb)/dy cross section due to errors of nPDFs and compares it with the Run 1 (upper panel) and Run 2 (lower panel) LHC data.The calculations using the central sets of nPDFs are given by the blue solid (EPPS21), red dashed (nCTEQ15WZSIH), and green dotted (nNNPDF3.0)curves and the error bands are represented by the corresponding shaded regions.One can see from the figure that within the PDF uncertainties the framework of NLO pQCD describes the data rather well; the agreement with the data is very good at central rapidity for all three nPDF sets (by construction), continues to be good for nCTEQ15WZSIH in the entire range of measured y, but becomes somewhat worse at higher |y| for EPPS21 and NNPDF3.0.
A comparison of our EPPS21 results with the previous EPP16 ones [24] shows that the full PDF uncertainty band, which receives contributions from varying the parameters of nPDFs and the baseline free proton PDFs, has come down to the order of few millibarns.As we discussed in Ref. [24], the free proton CT14nlo PDFs accompanying the EPPS16 nPDFs contain a particular error set dramatically growing at small x, which results in an abnormally large small-x uncertainty.In the new EPPS21 nPDFs, where the nuclear effects are correlated with the baseline CT18ANLO [5] free proton PDF error sets, this behaviour no longer persists.
One can clearly see from Fig. 4 that the EPPS21 nPDFs correspond to significantly smaller uncertainties than nNNPDF3.0 and CTEQ15WZSIH.In particular, the nNNPDF3.0uncertainties, which also account for the free proton PDF errors, at central rapidity rise up to around 5.6 mb at Run 1 and up to around 9.5 mb at Run 2. The nCTEQ15WZSIH uncertainties, which account for the nPDF errors only, are smaller both for Run 1 and Run 2 at central rapidities than at y ≈ ±2.0, which leads to a valley-like structure.For instance, for Run 2, the uncertainty rises up to around 8.5 mb at y = 0 and then to its maximum of approximately 18 mb at y ≈ ±2.5.As with EPPS21, no single error PDFs set stands out in the nCTEQ15WZSIH and nNNPDF3.0parametrizations, but the larger uncertainty bands are simply the result of a wider distribution in the underlying error sets. is approximately 1,000 times smaller than that in the Pb-Pb case primarily due to the much smaller photon flux.On the other hand, the shape of the y dependence is similar in the O-O and Pb-Pb cases: it is rather broad at midrapidity with sloping "shoulders" at forward and backward rapidities; higher scales correspond to larger dσ(O+O → O+J/ψ +O)/dy, which also tend to develop a valley-like structure at the highest scales of µ ≈ M J/ψ .
To quantify the magnitude of the scale dependence, we consider the ratio between the µ = M J/ψ and µ = m c results at y = 0 which we denote by R scale .One can see from Fig. 5 that R scale is of the same order of magnitude as in Pb-Pb collisions starting at R scale ≈ 16 at √ s N N = 2.76 TeV and rising up to R scale ≈ 35 at √ s N N = 7 TeV.We have checked that with nCTEQ15WZSIH the scale dependence is of the same order as with EPPS21: huge scale dependence is due to the nearly perfect cancellation between the LO and the NLO contributions in both the real and the imaginary parts of the amplitude at the lowest scale of µ = m c .At forward and backward rapidities over the full range µ ∈ [m c , M J/ψ ], the scale dependence is not as strong for all three nPDF sets under consideration.Figure 6 shows the separate contributions of the two terms to dσ(O + O → O + J/ψ + O)/dy in Eq. ( 1), labeled "W + " (dashed orange) and "W − " (dotted green), along with their sum labeled "Full" (solid blue).The calculation is carried out using the EPPS21 nPDFs at µ = 2.39 GeV, which is the optimal scale in the Pb-Pb case.The results are qualitatively similar to those for the Pb-Pb collision system [24].Looking only at the W + contribution, we observe a small bump at backward rapidities caused by the interplay of the large photon flux with the increasing photoproduction cross section and the integral of the nuclear form factor squared.This increase in the differential cross section is momentarily halted and then decreases as one moves from y ≈ −4 to y ≈ −2 (i.e. at Run 1 √ s N N and slightly differently for the other energies).Then the growth of the photopro-duction cross section forces an increase of the absolute magnitude of the UPC cross section until around y ≈ 2, when the decrease in the photon flux eventually forces the cross section to zero.One can see that this holds for all four energies and we have checked that the results are qualitatively similar for all the three nPDF sets studied here.
Figure 7 quantifies the contributions of the imaginary and real parts of the γ + A → J/ψ + A amplitude to the dσ(O + O → O + J/ψ + O)/dy UPC cross section: the dashed orange curve gives the result, when only the imaginary part is included, the dotted green curve shows the result, when only the real part is included, and the solid blue curve is their sum.One can see from the figure that with increasing √ s N N , the imaginary part becomes more important at central rapidity and, when moving from 2.76 TeV to 6.37 TeV, the dip in the imaginary part at around y ± 3 seen at √ s N N = 2.76 TeV actually rises above the real part, i.e., the imaginary part becomes the dominant contribution at all values of rapidity.Qualitatively, the results are the same for the other two nPDF sets nNNPDF3.0 and nCTEQ15WZSIH.Finally, in Fig. 8 we show the separate contributions of different parton channels to the UPC cross section.The dashed orange curve gives the gluon contribution, i.e., it corresponds to the situation when the contribution of quarks is neglected, the dotted green line gives the quark contribution, the red dash dotted curve is the interference  gluons begin to dominate at forward and backward rapidities.We have checked that this trend also persists for the nNNPDF3.0 and nCTEQ15WZSIH nPDFs.Lastly, a few words about the feasibility of measurements of this process in O-O UPCs.Experimentally the dσ coh J/ψ /dy rapidity differential cross section for the coherent photoproduction of J/ψ in the lepton channel l + l − is given by [38] where N coh J/ψ is the yield, i.e., the number of observed J/ψ particles, E is the combined acceptance and efficiency of the detector, Γ l + l − is the branching ratio to the desired final state l + l − , L int is the integrated luminosity, and ∆y is the width of the rapidity interval under consideration.By considering only the central rapidity and the the muon channel with Γ l + l − = 5.961% [47] and taking the values given in [38], E = 4.57 %, ∆y = 1.8, and N coh J/ψ = 250, together with dσ coh J/ψ /dy = 2 µb from Fig. 6, we can estimate the required integrated luminosity L int to be It was discussed in Ref. [25] that in the high luminosity O-O run at the LHC, the average luminosity would be ⟨L AA ⟩ = 8.99 × 10 30 cm −2 s −1 .This means that in a specialized 24-hour O-O run at ALICE, the integrated luminosity would be approximately 7.8×10 5 µb −1 resulting in approximately 7.5 × 10 3 J/ψ's making the experimental data acquisition more than feasible.Unfortunately, at the proposed short data acquisition during Run 3, one would most likely acquire only the integrated luminosity of 500 µb −1 , which means that one expects to see only five events [25].generates an additional scale uncertainty due to the fact that the O-O process will be probed at a smaller x value due to the skewness parameter ξ becoming smaller.
Figures 9, 10 and 11 present our NLO pQCD predictions for R O/Pb evaluated at six different values of the scale µ ranging from µ = 1.55 GeV to µ = 3.1 GeV using the EPPS21, nNNPDF3.0 and nCTEQ15WZSIH nPDFs, respectively.One can see from the figures that the relative scale uncertainty seems to be the smallest for EPPS21 and nCTEQ15WZSIH at y ≈ 0, which then grows slightly towards backward and forward rapidities.However, in the nNNPDF3.0case the situation is reversed due to the almost exact cancellation of the photoproduction amplitude for the O-O process at central rapidity.Moreover, depending on the energy, the EPPS21 nPDF set produces a node at y ≈ ±1.1 or y ≈ ±1.8, where all the scales except for the lowest µ = m c seem to agree with each other.Such a node is missing in the results given by nNNPDF3.0 or nCTEQ15WZSIH.In addition, we would like to point out that our predictions for R O/Pb for each nPDF set separately tend to cluster together at higher values of µ.
To quantify the magnitude of the relative scale dependence, we consider the super-ratio of ratios R O/Pb at y = 0, which are evaluated at µ = M J/ψ and µ = m c , The results for R O/Pb scale are presented in Table I.One can see from the table that for all three sets of nPDFs, the scale uncertainty of R O/Pb scale is smaller by approximately a factor of 10 than that of the predictions for the individual Pb-Pb and O-O UPC cross sections (the exact size of the reduction in the scale dependence depends on the particular nPDF set and √ s N N ).The scale uncertainty also increases, when √ s N N is increased, since at higher energies one probes the nPDFs at progressively smaller x, where the scale evolution of the nPDFs is faster.One can see from the table that the scale uncertainty characterized by the ratio R O/Pb scale of Eq. ( 17) turns out to be very large in the case of nNNPDF3.0 nPDFs.This √ s N N , while only moderately affecting the EPPS21 and nCTEQ15WZSIH results.
To better understand the scale and energy dependence of the ratio of the O-O and Pb-Pb UPC cross sections, we consider the ratio R O/Pb , when the numerator of Eq. ( 16) -the dσ(O + O → O + J/ψ + O)/dy cross section -is evaluated at √ s N N = 6.37 and 7 TeV, and the denominator of Eq. ( 16) -the dσ(Pb + Pb → Pb + J/ψ + Pb)/dy cross section -is evaluated at √ s N N = 2.76 TeV (Run 1) and 5.02 TeV (Run 2).Our results for the EPPS21, nNNPDF3.0 and nCTEQ15WZSIH nPDF sets are presented in Figs. 12, 13 and 14, respectively.One can see from the figures that qualitatively the scale dependence is similar for nNNPDF3.0 and nCTEQ15WZSIH, but in the case of EPPS21 for Pb-Pb UPCs at √ s N N = 2.76 TeV the node disappears and the systematics of the scale dependence becomes similar for all three nPDF sets.The general effect of taking R O/Pb scale at different energies means that the scale dependence is increased as given in Table II.As we take O-O consistently at a higher energy, dσ/dy increases scale by scale, as was shown in Fig. 5.For EPPS21 the situation is more involved since for the first two entries the scale dependence is flipped, but the magnitude of the dependence stays the same.For the last two entries -i.e., Pb-Pb taken at √ s N N = 5.02 TeV -the dependence actually gets smaller.For nNNPDF3.0 the situation is the worst: taking the ratio at different energies means that we increase the scale dependence by a factor of three at worst.For nCTEQ15WZSIH the fac-tor is only about 1.6.However, if we disregard the lowest scale and take µ = 1.74 GeV instead, the scale dependence becomes smaller for all three sets.For nNNPDF3.0 the drop is quite sizeable again: at all energies the scale dependence drops to less than a factor of 3. Figure 15 illustrates the PDF uncertainties of our NLO pQCD predictions for R O/Pb as a function of y for EPPS21, nNNPDF3.0 and nCTEQ15WZSIH nPDFs.The calculations using the central sets of the nPDFs at their corresponding optimal scales are given by the blue solid (EPPS21), green dotted (nNNPDF3.0),and red dashed (nCTEQ15WZSIH) curves.The corresponding uncertainties are given by the shaded bands.They are calculated by first finding the ratio R O/Pb for each error set and then using the asymmetric form (see Eq. ( 13)) for EPPS21 and nCTEQ15WZSIH and the CL prescription for nNNPDF3.0.Thus, the blue and green bands give the full (free proton + nuclear) uncertainty of EPPS21 and nNNPDF3.0nPDFs, respectively, while the hashed red band gives the nuclear uncertainty of the nCTEQ15WZSIH set.We see that the predictions agree within the PDF uncertainty bands.
One can see from the figure that as |y| is increased, the uncertainty bands grow bigger for all three sets.It can be understood by noticing that at higher positive  A comparison of the PDF and scale uncertainties in R O/Pb at y = 0 as a function of √ s N N is shown in Fig. 16.The PDF uncertainties are calculated at the corresponding optimal scales for the EPPS21 (left), nNNPDF3.0(middle) and nCTEQ15WZSIH (right) nuclear PDFs.The scale uncertainty represents the range between the scales µ = m c and µ = M J/ψ .In absolute terms the EPPS21 PDF uncertainty is typically smaller than the scale uncertainty, while for nNNPDF3.0 and nCTEQ15WZSIH the scale and PDF uncertainties are of similar magnitudes.The figure also shows the lack of uniformity between the uncertainties between different sets.For instance, in the EPPS21 case, the scale uncertainty dominates upwards, whereas the PDF uncertainties dominate downwards.For nNNPDF3.0, the situation is different: the scale uncertainties dominate the downwards uncertainty, while the PDF uncertainties are symmetric.Then interestingly for nCTEQ15WZSIH -the set with the enhanced strange quark contribution -the scale uncertainties are smaller than the PDF uncertainties at all energies.The value of the ratio stays approximately constant as a function of √ s N N for all three sets.Figure 17 presents the nPDF uncertainties of the ratio R O/Pb as a function of y, when the O-O and Pb-Pb UPC cross sections are evaluated at different collision energies (see our discussion above).The notation of the curves and shaded bands is the same as in Fig. 15.A comparison with Fig. 15 shows that the results in the two figures are similar.In particular, at central rapidity for EPPS21 the ratio between the upper bound and the lower bound for the PDF uncertainties is about 2.2 for Pb-Pb taken at Run 1 energy and 1.8 for Pb-Pb taken at Run 2 energy, which means that the PDF uncertainty is slightly larger at all energies under consideration than the scale uncertainty (see Table II).For nNNPDF3.0 the same ratio is around 3 for all energies and again the scale uncertainty is clearly the dominating one, when considering µ ∈ [m c , M J/ψ ].If we ignore the lowest scale µ = m c , we find that the PDF uncertainty is again the larger one.For nCTEQ15WZSIH the corresponding ratios are about 1.7 and 2.5 for Run 1 and Run 2 energies, respectively.
These results are summarized in Fig. 18, which shows the R O/Pb ratio at y = 0 for EPPS21, nNNPDF3.0, and nCTEQ15WZSIH for different configurations of collision energies as discussed above.The color-coded bars give the scale (wide error bars) and PDF (thin error bars with caps) uncertainties; the former are calculated using the central sets of the respective nPDF fits and the latter are

IV. CONCLUSIONS
This work continues our studies of J/ψ photoproduction in nucleus-nucleus UPCs at the LHC within the framework of collinear factorization and NLO perturbative QCD.In particular, we update our results for this process in Pb-Pb UPCs and make predictions for the dσ(Pb + Pb → Pb + J/ψ + Pb)/dy cross section as a function of the J/ψ rapidity y using the stateof-the-art EPPS21, nNNPDF3.0, and nCTEQ15WZSIH nPDF sets.Taking nuclear generalized parton distribution functions in their forward limit, where they reduce to the nPDFs, we obtain a good description of Run 1 and Run 2 LHC data on dσ(Pb + Pb → Pb + J/ψ + Pb)/dy.This is achieved by choosing an "optimal" scale for each set of nPDFs: µ = 2.39 GeV for EPPS21, µ = 2.22 GeV for nNNPDF3.0 and µ = 2.02 GeV for nCTEQ15WZSIH.
Compared to our earlier calculations using EPPS16, nNNPDF2.0, and nCTEQ15 nPDFs [24], we can make the following observations.The results employing the central set of the EPPS21 nPDFs are found to be similar to those with the EPPS16 nPDFs with the corresponding "optimal" scale µ = 2.37 GeV.In addition, with the EPPS21 set, the PDF uncertainties have reduced significantly.At the same time, our results with the nNNPDF3.0nPDFs exhibit a much more regular behavior than those corresponding to the nNNPDF2.0nPDFs and, as a result, better reproduce the data.This is due to the fact that the gluon distribution in nNNPDF3.0grows at small x much slower than that in nNNPDF2.0nPDFs.The best description of the data at both central and forward/backward rapidities at Run 1 and Run 2 energies is achieved with the nCTEQ15WZSIH, which also performs better than the nCTEQ15 set.This is due to the strongly enhanced strange quark content at small x.Thus, at least at NLO pQCD, this process is a potential probe of the elusive strange quark PDFs.It should be kept in mind that, because the scale dependence is significant, the situation may still change at NNLO.
We also made detailed predictions for the O + O → O + J/ψ + O rapidity-differential UPC cross section in anticipation of the planned oxygen run at the LHC.Comparing with Pb-Pb UPCs, we observe that the shape of the rapidity distribution in the O-O case is qualitatively similar to that in Pb-Pb, but the former begins to develop  a valley-like structure around y = 0 at high enough scales µ ∼ M J/ψ .At central rapidity, the scale dependence of our results for O-O corresponding to the EPPS21 and nCTEQ15WZSIH nPDFs is slightly smaller than that for Pb-Pb collisions, but it is still of the same order of magnitude.For nNNPDF3.0, the situation is worse: the scale uncertainty grows to the order of 10 3 due to the nearly perfect cancellation of the sum of the LO and the NLO contributions both in the real part and the imaginary parts of the amplitude at the smallest scale of µ = m c .The decomposition of the O-O results into the W ± components, the imaginary and the real parts, and the gluon and quark contributions did not differ significantly from the results in the Pb-Pb case [24].Namely, the W ± contributions exhibit a two-bump structure; the imaginary part gives the dominant contribution over a larger range of y, while the real part cannot be neglected, especially for large values of |y|; the quark contribution dominates at central rapidity, but the gluons become important at backwards or forwards rapidities.Furthermore, the interplay between the gluon and the quark contributions plays an important role.
In order to reduce the significant scale and nPDF uncertainties, we have studied the ratio of the J/ψ rapidity distributions in O-O and Pb-Pb UPCs at different collision energies √ s N N .We found that for EPPS21 and nCTEQ15WZSIH, the scale uncertainties in the ra-tio indeed became significantly smaller.The reduction in the scale dependence is largest at central rapidities and slightly smaller towards backward and forward rapidities both when the ratio is taken at the same value of √ s N N and when taken at different values.For nNNPDF3.0, the situation is the same at central rapidity, i.e., the ratio has a smaller scale dependence when compared to the O-O case.Interestingly, and contrary to EPPS21 and nCTEQ15WZSIH, the scale dependence of the nNNPDF3.0ratio at forward and backward rapidities becomes even smaller since the LO and NLO contributions in the O-O results no longer cancel to such an exact degree.
The PDF uncertainties for the ratios of the rapidity distribution in O-O to Pb-Pb UPCs for EPPS21 were found to be smaller than those for nCTEQ15WZSIH and nNNPDF3.0.This is a direct consequence of the tightly constrained error sets in EPPS21, whereas in nCTEQ15WZSIH and nNNPDF3.0,there is more variation.The comparison of the PDF and scale uncertainties for the ratios taken at the same energy shows that the scale uncertainty is the dominant one for the EPPS21, while for nCTEQ15WZSIH the situation is reversed.For nNNPDF3.0, the two uncertainties are similar.For the ratios taken at different energies, the PDF uncertainties are of the same magnitude for EPPS21 and nCTEQ15WZSIH and somewhat larger for nNNPDF3.0.Our analysis demonstrates that the large scale uncertainty of our NLO pQCD results can be tamed through suitably considered ratios of rapidity differential cross sections.In future work, it would be instructive to extend our analysis to J/ψ photoproduction in p-Pb and p-O asymmetric UPCs and also to photoproduction of Υ mesons in nucleus-nucleus UPCs.In addition, our framework could be improved through a more detailed GPD modeling [48] and the inclusion of the non-relativistic QCD corrections to the charmonium wave function [49][50][51][52].

FIG. 1 .
FIG. 1.The scaled photon flux (1/Z 2 )kdN A γ (k)/dk as a function of the rapidity y in the plus direction for the oxygen and lead beams for four different values of the c.m.s.energy √ sNN = 2.76, 5.02, 6.37 and 7 TeV.

CFIG. 4 .
FIG.4.The PDF uncertainties of the NLO pQCD predictions for the dσ(Pb + Pb → Pb + J/ψ + Pb)/dy cross section as a function of y for Run 1 (upper) and Run 2 (lower) at the LHC, and a comparison with the corresponding Run 1[38][39][40] and Run 2[41][42][43][44] data, mirrored with respect to y = 0 and with the statistical and systematic errors added in quadrature.The results corresponding to the central sets of nPDFs are shown by the blue solid (EPPS21), red dashed (nCTEQ15WZSIH), and green dotted (nNNPDF3.0)curves, respectively, and the error bands are represented by the corresponding shaded regions.All calculations are performed at the indicated values of the optimal scale µ.

FIG. 5 .
FIG. 5.The NLO pQCD results for the rapidity differential cross section of coherent J/ψ photoproduction in O-O UPCs as a function of the rapidity y, obtained with the EPPS21 nPDFs at √ sNN = 2.76, 5.02, 6.37 and 7 TeV.The different lines show the results for ten choices of the scale µ ranging from µ = mc (lowest curve) to µ = M J/ψ (highest curve) with a step of mc/8.The µ = 2.39 GeV "optimal scale" prediction lies in the middle of this scale-uncertainty envelope.

FIG. 6 .
FIG. 6. Separation of the NLO pQCD predictions for the dσ(O + O → O + J/ψ + O)/dy cross section of coherent J/ψ photoproduction in O-O UPCs as a function of the rapidity y into the W + (dashed orange curve) and W − (dotted green curve) components; the solid blue line is their sum.The calculation employs the EPPS21 nPDFs at µ = 2.39 GeV.The different panels correspond to √ sNN = 2.76, 5.02, 6.37 and 7 TeV.

FIG. 7 .
FIG. 7. The contributions of the imaginary (dashed orange curve) and real (green dotted curve) parts of the γ + A → J/ψ + A amplitude to the dσ(O + O → O + J/ψ + O)/dy cross section of coherent J/ψ photoproduction in O-O UPCs as a function of the rapidity y; the solid blue curve is the full result.The calculation uses the EPPS21 nPDFs at µ = 2.39 GeV.The different panels correspond to √ sNN = 2.76, 5.02, 6.37 and 7 TeV.

D.
Ratios of O-O and Pb-Pb UPC cross sectionsOur results presented above indicate that the scale dependence is considerable for both O-O and Pb-Pb collision systems.To reduce it, we examine the following scaled ratio of the O-O and Pb-Pb UPC cross section,R O/Pb = 208Z Pb 16Z O 2 dσ(O + O → O + J/ψ + O)/dydσ(Pb + Pb → Pb + J/ψ + Pb)/dy (16) where the factor of [(208Z Pb /(16Z O )] 2 is introduced to remove the effects of the Z 2 scaling of the photon flux and the A 2 scaling of the nuclear form factor squared.Since the hard scattering part is the same for both O-O and Pb-Pb scatterings, the scale dependence, which we expect to see in this ratio, comes from the underlying

FIG. 9 .
FIG. 9.The NLO pQCD predictions using the EPPS21 nPDFs for the scaled ratio of cross sections of J/ψ photoproduction in O-O and Pb-Pb UPCs as a function of the rapidity y for six different values of the scale µ at four different values of √ sNN .
FIG. 12.The scaled ratio of the NLO pQCD cross sections of J/ψ photoproduction in O-O and Pb-Pb UPCs as a function of the rapidity y for six different values of the scale µ at non-equal values of O-O and Pb-Pb collision energies.The results are obtained with the EPPS21 nPDFs.

FIG. 15 .
FIG. 15.The PDF uncertainties of NLO pQCD predictions for R O/Pb as a function of the rapidity y.The results corresponding to the central nPDF sets at the optimal scales are shown by the blue solid (EPPS21), green dotted (nNNPDF3.0),and red dashed (nCTEQ15WZSIH) curves, respectively.The corresponding uncertainties are shown by the shaded bands, see text for details.Different panels correspond to different √ sNN .

FIG. 16 .
FIG. 16.Comparison of the PDF (thin blue) and scale (wide orange) uncertainties in the ratio of the NLO pQCD calculation of O-O to Pb-Pb rapidity differential cross section at central rapidity, y = 0, for three different nPDF sets: EPPS21 (left), nNNPDF3.0(center) and nCTEQ15WZSIH (right).Here O-O and Pb-Pb are taken at same energy and all sets at their corresponding optimal scales.
FIG. 17.The scaled ratios for EPPS21 (solid blue), nNNPDF3.0(dotted green) and nCTEQ15WZSIH (dashed red) at their optimal scales as a function of the J/ψ rapidity.The blue band gives the EPPS21 uncertainty, the green band gives the nNNPDF3.090% CL uncertainty and the hatched red band gives the nCTEQ15WZSIH nuclear uncertainty.In the first row Pb-Pb has been taken at Run 1 energy and in the second row at Run 2 energy.The O-O energies correspond to the two proposed energies of 6.37 TeV (left column) and 7 TeV (right column).

FIG. 18 .
FIG. 18.The scaled ratios of O-O to Pb-Pb rapidity differential cross sections for EPPS21, nNNPDF3.0 and nCTEQ15WZSIH at their corresponding optimal scales at central rapidity, y = 0, where O-O and Pb-Pb have been taken at different √ sNN energies.In the left panel the Pb-Pb collision is taken at Run 1 energy and in the right panel at Run 2 energy.