Nuclear suppression from coherent $J/\psi$ photoproduction at the Large Hadron Collider

Using the data on coherent $J/\psi$ photoproduction in Pb-Pb ultraperipheral collisions (UPCs) obtained in Runs 1 and 2 at the Large Hadron Collider (LHC), we determined with a good accuracy the nuclear suppression factor of $S_{Pb}(x)$ in a wide range of the momentum fraction $x$, $10^{-5} \leq x \leq 0.04$. In the small-$x$ region $x<10^{-3}$, our $\chi^2$ fit favors a weakly decreasing or flat form of $S_{Pb}(x) \approx 0.6$ with a $10-15$% error at $x=5 \times 10^{-4}$ and a $20-30$% error at $x=10^{-4}$. At large $x$, $S_{Pb}(x)$ is constrained to 10% accuracy up to $x=0.04$ and is also consistent at $\langle x \rangle =0.042$ with the Fermilab data on the $A$ dependence of the cross section of coherent $J/\psi$ photoproduction on fixed nuclear targets. The resulting uncertainties on $S_{Pb}(x)$ are small, which demonstrates the potential of the LHC data on coherent charmonium photoproduction in Pb-Pb UPCs to provide additional constraints on small-$x$ nPDFs.


Introduction
Determination of nuclear parton distribution functions (nPDFs) is an important topic of phenomenology of high energy nuclear physics. In the context of Quantum Chromodynamics (QCD), collinear nPDFs are universal quantities encoding the microscopic quark and gluon structure of nuclei probed in various hard processes. One usually determines nPDFs using so-called global QCD fits to available data [1,2,3,4,5,6,7,8]. However, because of limited kinematic coverage of the available data and largely indirect determination of the gluon distribution, nPDFs are currently known with large uncertainties. Recent QCD analyses of the proton-nucleus (pA) data collected during Runs 1 and 2 at the LHC [9,10,11,12,13,14,15,16] showed that while they provide certain new restrictions on nPDFs, the remaining uncertainties are still significant. Alternatively, nuclear structure functions and nPDFs at small x can be theoretically predicted using models of nuclear shadowing, which are based on its connection to diffraction [17,18,19,20]. In particular, in Ref. [20], the use of QCD factorization theorems allowed one to connect the leading-twist nuclear shadowing of nPDFs to proton diffractive PDFs and, hence, predict small-x nPDFs with a small uncertainty. This topic will be further pursued after luminosity and energy upgrades of the LHC [21].
In the limit of very high energies, one often uses the framework of the color dipole model, which allows one to study the proximity of the dipole-target interaction to the new QCD regime characterized by saturation of the gluon density, for reviews, see, e.g. [22,23]. Establishing the pattern and signs of the saturation using HERA and LHC data remains a challenge and an active field of research, see, e.g. [24,25].
In the future, it is expected that nPDFs and possible signs of an onset of the gluon saturation will be explored with high precision and in a broad kinematic range using such lepton-nucleus colliders as the Electron-Ion Collider in the USA [26,27] and the Large Hadron-Electron Collider (LHeC) [28] and Future Circular Collider (FCC) [29] at CERN. Meanwhile it is important to utilize all existing capabilities of the LHC to constrain nPDFs including those provided by ultraperipheral collisions (UPCs) of heavy ions.
In UPCs, ions in colliding beams interact at large distances between their centers in the transverse plane (large impact parameters) so that strong hadron interactions are suppressed leading to the dominance of long-distance electromagnetic processes induced by ultrarelativistic nuclei, which in the equivalent photon approximation are characterized by fluxes of quasireal photons of high intensity and energy. Thus, it gives an opportunity to study photon-nucleus scattering and nPDFs at unprecedentedly high energies [30]. In particular, QCD analyses of photoproduction of heavy quarkonia [31,32,33] and inclusive and diffractive dijet photoproduction [34,35,36] at the LHC provided new information on nuclear gluon and quark distributions at small x.
This work continues and extends our phenomenological studies of nuclear suppression in coherent J/ψ photoproduction on nuclei at the LHC [32,33] by including in the analysis all the data available to date on the rapidity y dependence of the cross section of coherent J/ψ photoproduction in Pb-Pb UPCs at √ s N N = 2.76 TeV [37,38,39] and √ s N N = 5.02 TeV [40,41]. As a cross-check and reference point at lower energies, we test our results against the data on the mass number A dependence of the cross section of coherent J/ψ photoproduction on fixed nuclear targets (Be, Fe, and Pb) obtained at Fermilab [42]. Note that it would also be very beneficial to collect high statistics on J/ψ photoproduction in heavy-ion UPCs at RHIC because it would cover the x range of x ∼ 0.015 at y = 0 and help to further constrain the results of our analysis, provided the data is enough accurate. Expressing our results in terms of the nuclear suppression factor of S P b (x), we determine S P b (x) with a good accuracy in a wide range of x, 10 −5 ≤ x ≤ 0.04. The resulting uncertainties are much smaller than those of nPDFs for these values of x, which demonstrates the potential of the LHC data on coherent charmonium photoproduction in Pb-Pb UPCs to provide new constraints on small-x nPDFs and possible signs of saturation of the gluon density.

Nuclear suppression factor for coherent J/ψ photoproduction on nuclei
In UPCs, both colliding nuclei serve as a source of quasi-real photons and a target. Therefore, using the method of equivalent photons [43,44], the cross section of coherent J/ψ photoproduction in symmetric Pb-Pb UPCs is given by a sum of the following two terms where y is the rapidity of J/ψ, N γ/A (W γp ) is the photon flux, and σ γA→J/ψA (W γp ) is the photoproduction cross section containing all details of the strong photon-nucleus interaction and production of J/ψ. Note that interference of the two terms in Eq. (1) is sizable only at very small values of the J/ψ transverse momentum [45] and hence can be safely neglected.
In the laboratory frame (coinciding with centre-of-mass system in our kinematics), the measured rapidity of J/ψ can be related to the invariant photon-nucleon energy W γp , where E A is the nuclear beam energy and M J/ψ is the mass of J/ψ. The ambiguity in W γp for y = 0 is a reflection of the presence of two terms in Eq. (1), where the first term corresponds to the right-moving photon source and the plus sign in Eq. (2) and the second term corresponds to the left-moving photon source and the minus sign in Eq. (2) (provided that y is defined with respect to the right-moving nucleus emitting the photon).
To avoid inelastic strong ion-ion interaction destroying the coherence condition, the photon flux in Eq. (1) is calculated as convolution over the impact parameter b of the flux of quasireal photons emitted by an ultrarelativistic charged ion N γ/A (ω, b) [43,44] with the probability not to have inelastic strong ion-ion interactions Γ AA where ω = W 2 γp /(4E A ) is the photon energy; σ N N is the total nucleon-nucleon cross section; is the so-called nuclear optical density, which is calculated using the Woods-Saxon (two-parameter Fermi model) parametrization of the nuclear density ρ A [46]. One should emphasize that the precise determination of the photon flux using Eq. (3) in a wide range of ω is essential for the analysis of the present work. The validity of the equivalent photon approximation and a model [47,48] generalizing Eq. (3) were successfully tested in electromagnetic dissociation with neutron emission in Pb-Pb UPCs [49].
The UPC cross section (1) is subject to nuclear modifications, which originate from the photon flux and the photoproduction cross section and which in general depend on the rapidity y and the collision energy √ s N N . To quantify the magnitude of nuclear corrections due to the strong dynamics encoded in the photoproduction cross section and to separate the two contributions in Eq. (1), it is convenient to introduce the nuclear suppression factor of S P b (x) by the following relation, see Refs. [32,33]: where The denominator in Eq. (4) is the coherent J/ψ photoproduction cross section in the impulse approximation (IA), where F A (t) is the nuclear elastic form factor and |t min | = x 2 m 2 N is the minimal momentum transfer squared (m N is the nucleon mass). In our work, F A (t) was calculated using the Woods-Saxon parametrization of the nuclear density [46]. The differential cross section of J/ψ photoproduction on the proton was parametrized in the form [32], which provides a good description of the available data at fixed targets [50,51,52] and at HERA [53,54], where C 0 = 342 ± 8 nb/GeV 2 , δ = 0.40 ± 0.01, W 0 = 100 GeV. For W γp ≤ 1 TeV, this parametrization is consistent with a power-law fit to the W dependence of the γp → J/ψp cross section extracted from the LHCb data on coherent J/ψ photoproduction in proton-proton UPCs at √ s N N = 7 TeV [55] and √ s N N = 13 TeV [56]. For higher photon energies W γp > 1 TeV, the extracted cross section shows a deviation from a pure power-law extrapolation of the HERA data, see the discussion in Ref. [56]. However, this region of W γp is not probed in the Pb-Pb UPCs data and, hence, does not affect the results of our analysis. Thus, the σ IA γA→J/ψA (W γp ) cross section is evaluated model-independently using data-driven parameterizations of the nuclear form factor and the γp → J/ψp differential cross section.
Introducing the UPC cross section in the impulse approximation dσ IA AA→J/ψAA /dy, one can present the square root of the ratio of the UPCs cross sections entering Eqs. (1) and (7) in the following form Without loss of generality, we will use y ≥ 0 and, hence, W + γp ≥ W − γp and x + ≤ x − . The advantage of Eq. (8) is that it relates the experimentally measured UPC cross section ratio on the left-hand side to the nuclear suppression factor of S P b (x) on the right-hand side. However, it involves S 2 P b (x) at two different values of x and is generally dominated by the , which complicates the separation of the x + and x − contributions and reliable extraction of the x + term corresponding to higher energies. Nevertheless, the use of all the available data on Pb-Pb UPCs collected during Runs 1 and 2 at the LHC along with a general parametrization of S P b (x) allows us to extract S P b (x) down to x ≈ 10 −5 with a good precision. Note that the two contributions to the UPC cross section can also be separated by measuring ion-ion UPCs accompanied by mutual electromagnetic excitation of colliding ions followed by forward neutron emission [57]. Unfortunately, the statistics of such measurements is currently too low.  Table 1 showing the rapidity intervals, the corresponding UPC cross sections dσ AA→J/ψAA /dy, and the values of x + and x − . The last column gives the [dσ AA→J/ψAA /dy)/(dσ IA AA→J/ψAA /dy)] 1/2 cross section ratio calculated using Eq. (8); the error is the sum of experimental statistical and systematic uncertainties as well as the 5% theoretical error on the IA cross section [32] added in quadrature.
To constrain the nuclear suppression factor of S P b (x) in as a broad range of x as possible, we apply Eq. (8) to the available UPC data listed in Table 1. Note that as explained in Refs. [32,33], the Run 1 ALICE data point at y = 0 [37] (the first entry in Table 1) unambiguously and modelindependently corresponds to The shape of S P b (x) as a function of x is unconstrained. In this work, we test the following simple piece-wise parametrization of S P b (x) where x 0 = 0.00112 and c ≥ 0. Our fit function contains three free parameters: a is determined by Eq. (9) and b and c are constrained by the low-energy W − γp and the high-energy W + γp contributions to the UPC cross sections, respectively. Note that we require that c ≥ 0, which is part of our model assuming that the factor of nuclear suppression is a monotonic function of x on the considered interval of 10 −5 < x < 0.05. The resulting S P b (x) as a function of x is illustrated in Fig. 1, where we show the fit to the Run 1 and LHCb Run 2 data (upper panel) and the fit to all the data in Table 1 (lower panel). The shaded bands represent the uncertainties of the respective fits. Also, we show the values of S P b (x = 0.00112) extracted from the Run 1 ALICE data at y = 0 [filled circle with the corresponding error, see Eq. (9)] and S P b (x = 0.042) determined using the fixed-target Fermilab data [open square with the corresponding uncertainty, see Eq. (12) and the discussion below].
Several features of the obtained results are noteworthy. First, all the fits favor either weakly decreasing or flat S P b (x) for x ≤ 0.001. It agrees with the small-x behavior of the g A (x, µ 2 )/[Ag p (x, µ 2 )] ratio of the nuclear and proton gluon distributions assumed in the EPS09 [3] and EPPS16 [7] nPDFs and is consistent within uncertainties with predictions of the leading twist model of nuclear shadowing [20]; see also the discussion in Sec. 3. Note that while the fit in the upper panel of Fig. 1 favors S P b (x) slightly decreasing with a decrease of x for x ≤ 0.001, the lower panel fit corresponds to a flat S P b (x) for x ≤ 0.001. This is a consequence of the fact that the Run 2 LHCb data points lie systematically lower than the Run 2 ALICE points, which moreover have smaller experimental errors, see Table 1 and Fig. 2. Second, the use of the Runs 1 and 2 data allowed us to obtain S P b (x) with a reasonable accuracy in a wide range of x. For instance, S P b (x) is determined with 10% error up to x = 0.01 − 0.04, with 10 − 15% error down to x = 5 × 10 −4 and 20 − 30% error down to x = 10 −4 [when the range of percentage values is given, the larger (smaller) value refers to the upper (lower) panel]. Note that the use of the Run 2 ALICE data in the fit leads to smaller uncertainties, which, however, comes with a price of the noticeably higher value of χ 2 . Figure 2 demonstrates how well the calculation of the cross section of coherent J/ψ photoproduction in Pb-Pb UPCs using Eq. (1) with the nuclear suppression factor of S P b (x) describes the available Run 1 (upper panel) and Run 2 (lower panel) LHC data. It shows dσ AA→J/ψAA ( √ s N N , y)/dy as a function for |y|, where the solid lines and the shaded band correspond to S P b (x) and its uncertainty from the lower panel of Fig. 1. One can see from the figure that within the experimental and fit uncertainties, one obtains a good description of the data. An examination of the lower panel of Fig. 2 demonstrates that the Run 2 ALICE data points lie systematically higher than the LHCb points. It results in a significantly higher value of χ 2 of the fit in the lower panel of Fig. 1 than that in the upper panel of Fig. 1.
In addition to the UPC data in Table 1, S P b (x) at large x can be further constrained using the Fermilab data on coherent J/ψ photoproduction by a photon beam with the average energy of 120 GeV on fixed nuclear targets of beryllium (Be), iron (Fe), and lead (Pb) [42], where the corresponding average value of x is x = 0.042. The measured yields are normalized to the incoherent cross section on Be and their nuclear mass number A dependence is fitted to the power law A α with α = 1.40 ± 0.06 ± 0.04. It is close to the expectation of the impulse approximation: a χ 2 fit to the A dependence given by Eq. (5) gives σ IA γA→J/ψA ∝ A 1. 44 . This indicates that nuclear corrections at these values of W γp and x are small.
To convert this result into the value of S P b (x) at x = 0.042, we calculate it using the optical limit of the Glauber model, see, e.g., Ref. [32], where σ V N (W γp ) is the charmonium-nucleon cross section, which we keep as a free parameter. While the application of such an approach to charmonium photoproduction at high energies is questionable, it provides an adequate estimate of σ J/ψN at the considered medium energy. It is known very well that the effect of nuclear shadowing encoded in Eq. (11) slows down the A dependence of hadron-nucleus cross sections: the stronger the nuclear absorption (the larger the value of σ V N ), the slower the A dependence of S A (x). In our analysis, we varied the value of σ V N so that the A dependence of the product σ IA γA→J/ψA S 2 A (x) reproduces that of the Fermilab data and found that σ V N ( W γp = 16.4 GeV) = 3 ± 3 mb. Note that this value agrees with the charmonium-nucleon cross section σ J/ψN obtained in the generalized vector dominance model and the coupled-channel generalized Glauber model framework [58] and in the QCD dipole formalism [59]. Substituting it in Eq. (11), we find that S P b ( x = 0.042) = 0.90 ± 0.10 .
As can be seen in Fig. 1, this value is consistent with the results of the fits to the LHC UPC data.
It is important to mention the analysis of [60], where the nuclear suppression factor of S P b (x) was extracted from measurements of coherent J/ψ photoproduction in ultraperipheral and peripheral Pb-Pb collisions at the LHC at  [60] reproduced the finding of [32]. At the same time, the dipole model generally predicts a somewhat smaller nuclear suppression [61,62,63,64], whose magnitude significantly depends on details of the model implementation including the choice of the charmonium wave function and the dipole cross section.

Implications for nuclear PDFs
As discussed in the Introduction, coherent J/ψ photoproduction in Pb-Pb UPCs at the LHC can be used to obtain new constraints on the nuclear gluon distribution at small x. Indeed, at the leading logarithmic approximation of perturbative pQCD and in the static limit for the charmonium wave function [65], there is direct correspondence between the suppression factor of S P b (x) and the ratio of the nuclear and nucleon gluon distributions R g (x, where κ A/N ≈ 0.9 − 0.95 is a small correction taking into account the slightly different dependence of the nuclear and proton gluon distribution on x. While the relation of Eq. (13) is subject to nextto-leading order (NLO) QCD corrections [66,67], a model-dependent relation between generalized parton distributions (GPDs) and usual parton distributions, and relativistic corrections to the charmonium wave function [68,69,70,71,72], it is nevertheless instructive to directly compare these two quantities. This is presented in Fig. 3, which shows S P b (x) and R g (x, µ 2 ) as functions of x. For the latter, we used the EPPS16 [7] and nCTEQ15 [5] nPDFs, and predictions of the leading twist model of nuclear shadowing [20]; all were evaluated at µ 2 = 3 GeV 2 [33]. Several features of the presented results deserve to be pointed out. First, in the entire studied range of x, 10 −5 < x < 0.005, the shapes and the magnitudes of S P b (x) and the EPPS16 R g (x, µ 2 ) are similar. At the same time, the uncertainty on R g (x, µ 2 ) is much larger than that on S P b (x), which indicates that the UPC data on coherent S P b (x) can potentially significantly reduce the current large uncertainty in the nuclear gluon distribution. Second, while the shapes of the EPPS16 and nCTEQ15 R g (x, µ 2 ) are similar, the latter corresponds to the stronger nuclear shadowing (suppression). Third, while predictions for the leading twist model agree with S P b (x) for x ≥ 10 −3 , the former predicts R g (x, µ 2 ), which noticeably decreases as x is decreased.

Conclusions
In this work, we analyzed the Runs 1 and 2 LHC data on coherent J/ψ photoproduction in Pb-Pb UPCs in terms of the nuclear suppression factor of S P b (x) using its generic parametrization and a χ 2 fit to the ratios of the measured UPCs cross sections to those calculated in the impulse approximation. It allowed us to determine S P b (x) with a reasonable accuracy in a wide range of x, 10 −5 ≤ x ≤ 0.04. In particular, in the small-x region x < 10 −3 , the fit favors a weakly decreasing or flat form of S P b (x) ≈ 0.6 with a 10−15% error at x = 5×10 −4 and a 20−30% error at x = 10 −4 . At large x, S P b (x) is constrained to 10% accuracy up to x = 0.04 and is also consistent with the value S P b (x = 0.042) = 0.90 ± 0.10, which we found from the A dependence of the cross section of coherent J/ψ photoproduction on fixed nuclear targets measured at Fermilab. The uncertainties in S P b (x) are small, which demonstrates the potential of the LHC data on coherent charmonium photoproduction in Pb-Pb UPCs to provide new constraints on small-x nPDFs and possible signs of saturation. It will also be very beneficial to collect high statistics on J/ψ photoproduction in heavy-ion UPCs at RHIC, which would cover x ∼ 0.015 at y = 0 and help to further constrain the results of our analysis, provided that the experimental accuracy is sufficiently high.