Study of exclusive photoproduction of charmonium in ultra-peripheral lead-lead collisions

The cross-sections of exclusive (coherent) photoproduction $J/\psi$ and $\psi(\mathrm{2S})$ mesons in ultra-peripheral PbPb collisions at a nucleon-nucleon centre-of-mass energy of $5.02\,\mathrm{TeV}$ are measured using a data sample corresponding to an integrated luminosity of $228\pm10\,\mathrm{\mu b}^{-1}$, collected by the LHCb experiment in 2018. The differential cross-sections are measured separately as a function of transverse momentum and rapidity in the nucleus-nucleus centre-of-mass frame for $J/\psi$ and $\psi(\mathrm{2S})$ mesons. The integrated cross-sections are measured to be \mbox{$\sigma^\mathrm{coh}_{J/\psi} = 5.965 \pm 0.059 \pm 0.232 \pm 0.262\,\mathrm{mb}$} and \mbox{$\sigma^\mathrm{coh}_{\psi(\mathrm{2S})} = 0.923 \pm 0.086 \pm 0.028 \pm 0.040\,\mathrm{mb}$}, where the first listed uncertainty is statistical, the second systematic and the third due to the luminosity determination. The cross-section ratio is measured to be \mbox{$\sigma^\mathrm{coh}_{\psi(\mathrm{2S})}/\sigma^\mathrm{coh}_{J/\psi} = 0.155 \pm 0.014 \pm 0.003$}, where the first uncertainty is statistical and the second is systematic. These results are compatible with theoretical predictions.


Introduction
Ultra-peripheral collisions (UPCs) occur when two nuclei collide with an impact parameter, the distance between their centres, larger than the sum of their radii [1]. Because the nuclei do not overlap, strong interactions are suppressed so that photon-induced interactions between the two ions dominate. The number of photons produced is proportional to the square of electric charge, so photon-nuclear interactions are significantly enhanced in leadlead (PbPb) collisions compared to proton-proton (pp) collisions. In UPCs, J/ψ and ψ(2S) mesons can be produced from the colourless exchange of a photon from one of the two nuclei and a pomeron from the other. Coherent (exclusive) photoproduction occurs when the photon couples coherently with the entire nucleus through an exchange of a pomeron, while for incoherent photoproduction, the photon interacts with a particular nucleon within the nucleus. In this work, the terms "coherent" and "incoherent" charmonium photoproduction refer to the two diagrams, respectively, shown in Figure 1. The coherent (exclusive) photoproduction of J/ψ and ψ(2S) mesons is expected to probe the nuclear gluon distribution functions at a momentum transfer of Q 2 ≈ m 2 /4, where m is the mass of the meson. The photon-nuclear production of these mesons depends on the longitudinal momentum fraction of gluons in the nucleus, x ≈ (m/ √ s NN )e ±y , where y is the rapidity of the meson and √ s NN is the nucleon-nucleon centre-of-mass energy.
Thus, coherent photoproduction of charmonium mesons provides an excellent laboratory to study nuclear shadowing effects and the initial states of collisions with small x, where 10 −5 ≲ x ≲ 10 −2 at the LHC [2]. The charmonia produced in this process have typical transverse momenta, p T , smaller than 100 MeV/c, with no other particles produced in the collision. Coherent (exclusive) J/ψ photoproduction was first measured in UPCs at HERA [3,4] in electron-proton scattering, and with ions at RHIC [5]. This process has also been measured by the CMS experiment in the central region [6], by the LHCb experiment in the forward region [7], and by the ALICE collaboration in both the central and forward regions [8,9] at √ s NN = 5.02 TeV in PbPb collisions at the LHC.
This paper presents a measurement of the coherent J/ψ and ψ(2S) production reconstructed through the dimuon final state using the 2018 PbPb data sample collected by the LHCb experiment at √ s NN = 5.02 TeV and corresponding to an integrated luminosity of 228 ± 10 µb −1 . The study also measures the ratio between the coherent ψ(2S) and J/ψ production cross-sections, where the uncertainties due to systematic effects and the luminosity determination largely cancel. This more precise measurement will help to constrain theoretical predictions, where uncertainties arise from the choice of the meson wave function in dipole scattering models [10,11] and the factorisation scale in perturbative QCD models [12].
The LHCb detector and simulation are described in Sec. 2. The selection of signal candidates and the determination of cross-sections are described in Sec. 3 and Sec. 4, respectively. The uncertainties due to systematic effects are described in Sec. 5, while the results are presented in Sec. 6 and conclusions are given in Sec. 7.

Detector, event reconstruction and simulation
The LHCb detector [13,14] is a single-arm forward spectrometer covering the pseudorapidity range 2 < η < 5, designed for the study of particles containing b or c quarks. The detector includes a high-precision tracking system consisting of a siliconstrip vertex detector surrounding the collision region, a large-area silicon-strip detector located upstream of a dipole magnet with a bending power of about 4 Tm, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet.
The tracking system provides a measurement of the momentum, p, of charged particles with a relative uncertainty that varies from 0.5% at low momentum to 1.0% at 200 GeV/c. The minimum distance of a track to a primary pp collision vertex, the impact parameter (IP), is measured with a resolution of (15 + 29/p T ) µm, where p T is in GeV/c. Photons, electrons and hadrons are identified by a calorimeter system consisting of scintillating-pad (SPD) and preshower detectors, an electromagnetic and a hadronic calorimeter. Muons are reconstructed as a long track passing through the vertex detector and the three stations of silicon-strip tracking detectors, and identified by a system composed of alternating layers of iron and multiwire proportional chambers.
The pseudorapidity coverage is extended by forward shower counters (HeRSCheL) consisting of five planes of scintillators with three planes at 114, 19.7 and 7.5 m upstream of the LHCb detector, and two planes downstream at 20 and 114 m. The HeRSCheL detector [15] significantly extends the acceptance for detecting particles from dissociated nucleons by covering the pseudorapidity range of 5 ≲ |η| ≲ 10, enhancing the classification of central exclusive production and UPC events.
The online event selection is performed by a trigger, which consists of a hardware stage, based on information from the calorimeter and muon systems, followed by a software stage, which applies a full event reconstruction.
Simulated events are used to determine corrections for the detector resolution, acceptance, and efficiency. The UPCs are modelled using STARlight [16] with a specific LHCb configuration [17]. The STARlight generator models coherent and incoherent vector-meson production in photon-nuclear interactions. Decays of unstable particles are described by EvtGen [18] with QED final-state radiation handled by Photos [19]. The interactions of the generated particles with the detector are modelled using the Geant4 toolkit [20,21] as described in Ref. [22].

Selection of signal candidates
Signal candidates are reconstructed through the decays J/ψ → µ + µ − and ψ(2S) → µ + µ − , and are required to have a rapidity within the range 2.0 < y * < 4.5, where the starred notation indicates that the observable is defined in the nucleus-nucleus centre-of-mass frame. All remaining selection criteria given here are defined in the laboratory frame. One of the candidate muons must pass the hardware-level trigger, which requires a muon p T greater than 500 MeV/c. The dimuon candidates are selected with a minimum-bias software trigger, requiring at least one track reconstructed by the vertex detector; this software trigger is 100% efficient with respect to the following offline selection, since it has a looser multiplicity requirement. The offline selection requires two muon candidates, both with tracks that have p T > 700 MeV/c within the pseudorapidity range 2.0 < η < 4.5. The dimuon candidates are required to have p T < 1 GeV/c and an azimuthal opening angle between the muons larger than 0.9 π. The mass of each signal candidate, m µ + µ − , is required to be within ±65 MeV/c 2 of the known J/ψ mass [23] or ±77.35 MeV/c 2 of the known ψ(2S) mass [23]. To suppress background from PbPb collisions with impact parameter smaller than two times the nucleus radius, only events with less than 20 hits in the SPD are retained, corresponding to very low occupancy events that make up about 0.3% of all minimum-bias events. Additionally, a requirement based upon a figure of merit that combines the signals from all HeRSCheL stations [15], is used to discard events with significant activity in the HeRSCheL acceptance region.

Cross-section determination
For comparison with theoretical predictions, the measured cross-sections are transformed into the nucleus-nucleus centre-of-mass frame, from the laboratory frame, to account for the non-zero crossing angle between the two Pb beams. The differential cross-section for coherent charmonium production in a given interval of rapidity or transverse momentum is determined as dσ coh where ψ is either J/ψ or ψ(2S), x represents either the y * or p * T , N coh ψ is the coherent signal yield, ε tot is the total efficiency, L is the integrated luminosity, ∆x is the width of either the y * -or p * T -interval, and B(ψ → µ + µ − ) is the branching fraction of the charmonium decay. The branching fractions B(J/ψ → µ + µ − ) = (5.961 ± 0.033) × 10 −2 and B(ψ(2S) → e + e − ) = (7.93 ± 0.17) × 10 −3 [23] are used. For the ψ(2S) the more accurate dielectron branching fraction is used, where lepton universality is assumed. The ratio between the differential cross-sections of ψ(2S) and J/ψ production in a given rapidity interval is given by The signal yields are extracted in two steps. First, an unbinned extended maximumlikelihood fit to the dimuon mass distribution is performed to obtain the J/ψ and ψ(2S) yields within the J/ψ and ψ(2S) mass windows, respectively. The nonresonant background yield, mostly due to γγ → µ + µ − process, is also obtained from the mass fit. This fit uses double-sided Crystal-Ball functions to describe the J/ψ and ψ(2S) mass shapes and an exponential function for the nonresonant background. The fit is performed in the range 2.9 < m µ + µ − < 4.0 GeV/c 2 . The mass distribution and the corresponding fit are shown in Fig. 2.
In the second step, the coherent yields are determined with unbinned maximumlikelihood fits to the ln(p * 2 T ) distributions separately for the candidates inside the J/ψ and ψ(2S) mass windows. The yields of J/ψ production include contributions from coherent and incoherent production, and feed-down from ψ(2S) decays into J/ψ (ψ(2S) → J/ψ X). Similarly, the ψ(2S) yields include contributions from both coherent and incoherent production, while the feed-down contribution from higher-order charmonium excited states is negligible given the current statistical precision. The quantity ln(p * 2 T ) is used because the variable p 2 T is a proxy for the typical momentum exchange, |t| ≈ p 2 T , in an elastic scattering process, and the logarithmic distribution allows one to see the peak of the data at low exchanged momenta. The coherent production has the smallest momentum exchange by definition, while the incoherent production gives a relatively larger transverse momentum to the J/ψ or ψ(2S) meson to balance the break-up of the pomeron-emitting nucleus. The feed-down contribution to J/ψ production also has greater transverse momentum than the coherent production to balance the other products from the ψ(2S) decay. The ln(p * 2 T ) shapes of coherent, incoherent and ψ(2S) feed-down components are taken from STARlight simulation, while the normalisation of these components are left free in the fit. The nonresonant background consists mostly of the γγ → µ + µ − process with a slightly lower transverse momentum of the dimuon system than coherent charmonium production. The distribution also contains a small contribution  from the random pairing of uncorrelated muons produced in the hadronic interactions during peripheral or central lead-lead collisions, signified by a large transverse momentum of the dimuon system. The STARlight simulation gives a precise description of ln(p * 2 T ) spectrum of the γγ → µ + µ − process, but not of the background from hadronic interactions. Instead, a data-driven method is chosen to model the nonresonant background by taking the dimuon candidates in the mass range 3.2 < m µ + µ − < 3.6 GeV/c 2 outside charmonium mass windows. In this way, the model includes the γγ → µ + µ − process and the QCD background together, and gives an unbiased modelling of the ln(p * 2 T ) spectrum. The yields of the nonresonant background are determined as the integral of the nonresonant component from the dimuon mass fit separately in the J/ψ and ψ(2S) mass windows, and are fixed in the ln(p * 2 T ) fits. Figure 3 shows the ln(p * 2 T ) distributions of selected J/ψ and ψ(2S) candidates in the rapidity interval 2 < y * < 4.5. Fits to the ln(p * 2 T ) distributions are performed in each y * interval to extract the corresponding J/ψ and ψ(2S) yields, as reported in Table 1. The coherent yield of J/ψ and ψ(2S) production for each p * T interval is calculated by subtracting the background components from the measured yield for that interval as reported in Tables 2 and 3. The contributions from background components are determined by an overall fit to the ln(p * 2 T ) distributions.  The total efficiency ε tot is determined as the product of the acceptance efficiency (ε acc ), the muon acceptance efficiency (ε µ-acc ), the tracking efficiency (ε trk ), the selection efficiency (ε sel ), the particle identification (PID) efficiency (ε PID ), the trigger efficiency (ε trg ) and the HeRSCheL-veto efficiency (ε her ). Each efficiency is evaluated separately for J/ψ and ψ(2S) mesons in each y * and p * T interval for the differential cross-section measurements. Efficiencies are evaluated from simulation calibrated to data. The value of ε µ-acc is determined at generator level as the fraction of events with both muon candidates passing p T > 700 MeV/c and 2.0 < η < 4.5. The signal candidates are required to pass the p T < 1 GeV/c selection and fall in the mass windows defined in Sec. 3, for J/ψ and ψ(2S) mesons separately. For ε trk , ε PID and ε trg , the simulation does not always describe the data well. Efficiency corrections from data using the tag-and-probe method [24] are determined from J/ψ events in PbPb collision data. The HeRSCheL-veto criteria is chosen to retain a signal efficiency of 90% according to a set of separately selected pure signal and background data samples. Dependencies of the efficiency correction on y * and p * T of the dimuon system are studied and found to be negligible in the invariant mass range from 2.9 to 4.0 GeV/c 2 .

Systematic uncertainties
Systematic uncertainties on the cross-section measurements arise from the efficiency and background determination, signal and background shapes, momentum resolution, integrated luminosity and knowledge of the J/ψ → µ + µ − and ψ(2S) → µ + µ − branching fractions. For the ψ(2S) to J/ψ cross-section ratio measurement, only systematic uncertainties from the charmonia decay branching fractions are considered. Those from efficiency and background determination, signal and background shapes integrated luminosity are highly correlated and cancel. A summary of the systematic uncertainties is presented in Table 4.
The systematic uncertainties related to the efficiencies are driven by the sizes of the simulation and data samples. They vary from (0.5 -2.0)% for the tracking efficiency, (0.9 -1.6)% for the PID efficiency and (2.1 -3.7)% for the trigger efficiency, across different y * and p * T intervals. The uncertainty associated with the HeRSCheL efficiency is a constant 1.4%.
The uncertainty on the background shape is estimated by varying the shape parameters within their fitted uncertainties. The maximum difference on the extracted signal yields is 1.2%, and is assigned as the background uncertainty.
The momentum resolution is expected to shift events from one p * T interval to another. The uncertainties due to the momentum resolution are evaluated by comparing the p T spectra between generated and reconstructed events. The evaluated relative uncertainties vary from 0.9 to 34% for different p * T intervals. The largest uncertainty corresponds to the p * T interval between 140 to 160 MeV as shown in Table 8 (Appendix A), where very small signal yields are observed.
The slight discrepancies between the data and the fit results are visible at ln(p * 2 T ) ∼ −4 [ln ( GeV/c 2 )] for both the J/ψ and ψ(2S) fits, as seen in Fig. 3. This is expected to originate from a mis-modelling of the predicted signal shape from simulation. A systematic uncertainty on the signal shape model is estimated by evaluating the difference between the fitted signal yields with respect to an alternative empirical signal shape. The obtained difference is about 0.04%, negligible compared to other uncertainties shown in Table 4.
The uncertainties on the branching fractions result in relative uncertainties on the measured cross-sections of 0.6% and 2.1% [23], respectively. The relative uncertainty on the luminosity is 4.4% [25].

Results and discussion
The integrated cross-sections of coherent J/ψ and ψ(2S) photoproduction in PbPb collisions are measured in the rapidity region 2.0 < y * < 4.5 as σ coh J/ψ = 5.965 ± 0.059 ± 0.232 ± 0.262 mb , σ coh ψ(2S) = 0.923 ± 0.086 ± 0.028 ± 0.040 mb , where the first listed uncertainty is statistical, the second is systematic and the third is due to the luminosity determination. The cross-section ratio between coherent ψ(2S) and J/ψ photoproduction is measured to be where the first uncertainty is statistical and the second is systematic. The luminosity uncertainty cancels in the ratio measurement. The measured differential cross-sections for coherent J/ψ and ψ(2S) photoproduction as functions of y * and p * T are shown in Figs. 4 and 5, respectively. The cross-section ratio of coherent photoproduction between ψ(2S) and J/ψ as a function of rapidity is shown in Fig. 6. The data are shown as black points with black error bars for the statistical uncertainties, red boxes show the systematic uncertainties and the fully correlated uncertainty due to integrated luminosity is labelled separately. In the same figures, the results are compared to several theoretical predictions. The numerical values of the results are reported in Tables 5 -9 in Appendix A.
The STARlight prediction is based on the concept of vector meson dominance with parameters tuned according to previous UPC data [16]. As shown in Figs. 4 and 5, it gives a good description of the decreasing slope as a function of y * and the shape as a function of p * T , but the overall predicted normalisation is about 20% and 50% higher for J/ψ and ψ(2S) production, respectively. The ratio between ψ(2S) and J/ψ production in Fig. 6 is also well modelled within data uncertainties.
Two sets of calculations using leading-order perturbative QCD (LO pQCD) are provided by Guzey, Kryshen, Strikman and Zhalov [12,26] (GKSZ) for both J/ψ and ψ(2S) coherent photoproduction. One uses the leading twist approximation (LTA) [27] to model the nuclear shadowing effect in the initial state. The shaded area labelled "LTA" in Fig. 4 corresponds to the uncertainties on the nuclear shadowing determined in Ref. [27]. The other uses EPS09 nuclear parton distribution functions (nPDFs) [28] for the nuclear shadowing, with an error band labelled "nPDF unce." under "EPS09" in Fig. 4 presenting the uncertainties of the nuclear modification. Note that the two LO pQCD calculations carry ad hoc normalisation factors of the cross-section determined using high-energy HERA data [12,29]. Both of them predict well the shapes of the data  for both J/ψ and ψ(2S) production as a function of y * in Fig. 4. A slightly larger (smaller) p * T is predicted for J/ψ (ψ(2S)) production than the data in Fig. 5. An underestimation of about 15% of the normalisation can be seen for both J/ψ and ψ(2S) production, but the ratio is well modelled in Fig. 6. The large nPDF uncertainties in Fig. 4 indicate that coherent charmonium photoproduction in heavy ion collisions is very sensitive to the nuclear modification factors, especially to the modelling of the gluon shadowing, used in the LO pQCD calculations [12].
The next-to-leading-order (NLO) pQCD calculation using the most recent EPPS21 NLO nPDFs [30] is provided by Flett, Eskola, Guzey, Löytäinen and Paukkunen [31] (FEGLP), and is only available for J/ψ production as shown in the left plot of Fig. 4. This is the first pQCD calculation without using ad hoc normalisation factors of the cross-section compared to previous LO calculations. The predicted central value is about 15 − 20% lower than the data, which is calculated based on a factorization/renormalisation scale, µ = 0.76 m J/ψ = 2.37 GeV, tuned using previous ALICE [8,9,32,33], LHCb [7] and CMS [6] data. The substantial shaded area labelled "scale variation" corresponds to a variation of µ from m J/ψ /2 to m J/ψ , indicating that the cross-section is extremely sensitive to the missing higher-order pQCD corrections. The nPDF uncertainties are much smaller for rapidity below 2 but much bigger for rapidity greater than 3 in the NLO pQCD calculation compared to the LO calculation. This is understood as an interplay of the real and imaginary parts of the quark and gluon amplitudes that causes a certain level of mutual cancellation of the nuclear effects, especially at lower rapidity region [31]. High-precision data can nevertheless help to further understand these effects.
Predictions calculated by Mäntysaari, Schenke and Lappi [38,39] (MSL) use the IP-SAT parameterisation to describe the dipole-proton cross-section but include sub-nucleon scale fluctuations. Calculations with (Is fluct.) and without (No fluct.) sub-nucleon fluctuation together with BG or GLC vector-meson wave functions are compared to the measurements in Figs. 4 and 5. Only the No fluct.+BG calculation is available for ψ(2S) as a function of y * and only No fluct.+BG and Is fluct.+BG are available for ψ(2S) as a function of p * T . They predict well the shapes of the J/ψ differential cross-section as functions of y * and p * T , but the predicted p * T for ψ(2S) production is slightly smaller than the data. Variations in the normalisation are relatively large among these models. The two models using BG vector-meson wave functions predict higher normalisation than the two using GLC. Because these models are calculated for rapidity below 3.5, the normalisation of the predictions as a function of p * T appears relatively lower than that as a function of y * . Calculations for ψ(2S) are relatively worse than for J/ψ because we know less about the ψ(2S) wave function than J/ψ [10,39,40], the precise data can nevertheless be helpful to improve this aspect. Among them, Is fluct.+BG gives the best prediction for J/ψ as a function of y * .
Models by Kopeliovich, Krelina, Nemchik and Potashnikova [35] (KKNP) are composed of quarkonium wave functions determined by the Buchmüller-Tye (BT) [43] or power-like (POW) [44,45] potentials, as well as the Golec-Biernat-Wusthoff (GBW) [46,47] or Kopeliovich-Schafer-Tarasov (KST) [48] models for the dipole-nucleon cross-sections. They appear similar to the models provided by GMMNS, with reasonably good agreement for J/ψ production for y * < 3, but with an overestimation of the decreasing slope, and consequently an underestimation of the data of 20-60% for y * > 3 for ψ(2S) production. But the decreasing slopes are consistent between the J/ψ and ψ(2S) calculations so that the predicted ratio between the two agrees well with the data.
In an alternative approach (GG-hs+BG) by Cepila, Contreras and Krelina [34] (CCK), the BG model is used for the vector-meson wave function, and the dipole-nucleon crosssection is parameterised assuming the nucleon is composed of so-called hot-spots (hs), regions with high-gluon density. The standard Glauber-Gribov (GG) formalism [49][50][51] is then used to extend the dipole-nucleon cross-section to the case for dipole-nucleus. This model describes well the slope as a function of y * for both J/ψ and ψ(2S) data, but a relatively large overestimation of the normalisation for ψ(2S) production. The corresponding prediction for the ratio between ψ(2S) and J/ψ production is therefore relatively higher than the data points.
In a closer look at the differential cross-section as a function of rapidity in Fig. 4 for both J/ψ and ψ(2S) mesons, one can observe that the data do not consistently decrease at a fixed slope, instead it has a subtle and elusive bump between 3 and 4. This is the first observation of this subtle signature thanks to the high precision data. Among the models discussed above, only the standard pQCD calculations can reproduce this feature, and can be understood as an interplay of the real and imaginary parts of the quark and gluon amplitudes [12,31].

Conclusion
The coherent (exclusive) J/ψ and ψ(2S) photoproduction cross-sections in PbPb ultraperipheral collisions at a centre-of-mass energy of √ s NN = 5.02 TeV are studied using a data sample corresponding to an integrated luminosity of 228 ± 10 µb −1 collected by the LHCb detector. The differential cross-sections, as a function of y * and p * T , are measured separately for J/ψ and ψ(2S) mesons in the ranges 2.0 < y * < 4.5 and 0 < p * T < 0.2 GeV/c. The ratio of the cross-sections between the coherent ψ(2S) and J/ψ production, as a function of rapidity, is also determined for the first time in PbPb collisions and is found to be compatible with theoretical models. The J/ψ results are the most precise measurement to date, while the ψ(2S) results represent the first measurement in the forward region.      Table 7: The differential cross-section ratio between ψ(2S) and J/ψ coherent production as a function of y * . The uncertainty due to luminosity determination cancels in the ratio.   LHCb collaboration