Probing subnucleon scale fluctuations in ultraperipheral heavy ion collisions

We show that introducing subnucleon scale fluctuations constrained by HERA diffractive $J/\Psi$ production data significantly affects the incoherent diffractive $J/\Psi$ production cross section in ultraperipheral heavy ion collisions. We find that the inclusion of the additional fluctuations increases the ratio of the incoherent to the coherent cross section approximately by a factor of $2$, and modifies the transverse momentum spectra of the produced $J/\Psi$ at momenta larger than the scale that corresponds to the distance scale of the subnucleonic fluctuations. We present predictions for $J/\Psi$ production in ultraperipheral heavy ion collisions at $\sqrt{s_{NN}}=5.02\,{\rm TeV}$ at the LHC and $200\,{\rm GeV}$ at RHIC.


I. INTRODUCTION
Deep inelastic scattering (DIS) is a clean process to study the internal structure of hadrons via interaction with (virtual) photons. The most precise data to date on the partonic structure of the proton comes from the DIS experiments performed at the HERA e ± + p collider [1,2]. These measurements have shown that the gluonic density of the proton grows rapidly when the momentum fraction x of the gluons decreases. Nonlinear QCD phenomena must limit this growth at very small x in order to avoid violating unitarity. These nonlinearities are most conveniently described within the Color Glass Condensate (CGC) effective theory of QCD [3,4], which defines the framework we use in this work.
The gluon density scales as A 1/3 , which amplifies nonlinear effects in heavy nuclei with large mass number A. However, currently there is no nuclear DIS data at small x available. The proposed Electron Ion Collider (EIC) [5,6] and Large Hadron Electron Collider (LHeC) [7] are designed to explore this region of large gluon densities. Before these machines are realized, one possibility to study deep inelastic scattering on nuclei is given by ultraperipheral heavy ion collisions, where a large impact parameter suppresses the relatively short range strong interactions. Instead, scattering processes involve a photon, emitted by one of the electrically charged nuclei, scattering off the other nucleus. First results from the LHC on diffractive vector meson production in ultraperipheral collisions have demonstrated the sensitivity of this process on nuclear effects at small x [8][9][10]. Recently, ultraperipheral proton-nucleus collisions, where the large electric charge of the nucleus causes the photon-proton scattering to dominate, have also been used to study deep inelastic scattering off a proton at energies much higher than were available at HERA [11,12].
Diffractive scattering, where a system of particles or just a single particle is produced without exchanging a net color charge with the target, is a powerful process to study the small-x structure of hadrons. At leading order in perturbative QCD, the diffractive cross section is proportional to the square of the gluon density, making it very sensitive to small-x gluons. In addition, exclusive vector meson production is sensitive to the geometric structure of the hadron. In particular, in coherent diffraction (where the target hadron remains in the same quantum state) an average density profile is probed. On the other hand, when the target breaks up, one is sensitive to the event-by-event fluctuations of the gluon fields in the target [13][14][15][16][17][18]. In recent phenomenological applications it has been demonstrated that one can indeed constrain the shape and the shape fluctuations of the proton (not included in previous literature) and nuclei at small x by studying exclusive vector meson production [19][20][21] (see also [22][23][24][25]). For a more formal introduction of diffractive scattering, the reader is referred to Ref. [26].
We calculate coherent and incoherent diffractive J/Ψ production in ultraperipheral heavy ion collisions in the dipole picture. The main focus of this work is to study how the diffractive cross sections are affected by inclusion of the subnucleon scale fluctuations that have been constrained using diffractive J/Ψ production at HERA in Refs. [20,32]. Apart from the inclusion of these fluctuations, we improve over previous state-of-the art CGC work [27] (see also [28][29][30][31]) by using a Monte Carlo method to explicitly calculate target averages. This allows us to use the impact parameter dependent saturation model (IPsat) dipole amplitude [33] without factorization of the impact parameter dependence (an approximation which was necessary to derive the incoherent cross section in Ref. [18]). We also use an updated IPsat parametrization fitted to the combined HERA data [34].

II. DIPOLE SCATTERING
The basic ingredient in the dipole framework is the forward elastic dipole-target scattering amplitude N (r, b, x P ), where r is the two dimensional vector that connects the quark and antiquark of the dipole in the transverse plane, b is the impact parameter and x P is the usual Bjorken variable of DIS in a diffractive event. In this work we use the IPsat model [33], which employs an eikonalized DGLAP evolved [35] gluon distribution xg and includes saturation effects. The dipole-proton scat-tering amplitude is written as (1) with the thickness function Here both the coupling α s and the gluon distribution are evaluated at the scale µ 2 = µ 2 0 + 4/r 2 . The proton width B p , initial scale µ 2 0 and the initial condition for the DGLAP evolution of the gluon distribution xg are parameters of the model. Their values are determined in Ref. [34] by performing fits to HERA DIS data.
To include proton structure fluctuations we follow Refs. [20,32] and assume that the gluonic density of the proton in the transverse plane is distributed around three constituent quarks (hot spots) 1 . These hot spots are assumed to be Gaussian. In practice we perform the replacement of the impact parameter profile (2) (3) where b q,i are the locations of the hot spots. They are sampled from a two dimensional Gaussian distribution whose width is B qc . The free parameters B q and B qc are obtained in Ref. [32] by comparing with the HERA coherent and incoherent diffractive J/Ψ production data at the photon-proton center of mass energy W ∼ 75 GeV, corresponding to x P ∼ 10 −3 . The proton fluctuation parameters obtained are B qc = 3.3 GeV −2 , B q = 0.7 GeV −2 and are close to the values obtained in a similar analysis in Ref. [21].
An additional source of fluctuations we include here comes from fluctuations of the overall normalization of the saturation scale, which we refer to in short as saturation scale fluctuations. Following again Ref. [32], we allow the saturation scale of each of the constituent quarks to fluctuate independently according to a log-normal distribution. The width of that distribution is obtained in Refs. [36,37] by comparing to the p+p multiplicity fluctuation data. The saturation scale fluctuations were shown in Ref. [32] to be necessary to describe the incoherent diffractive cross section measured by HERA at small |t|. For a more detailed description of the implementation see Ref. [32].
The replacement (3) changes the impact parameter dependence of the average dipole amplitude (1) even though Here, the average is taken over different nucleon configurations. This is because the thickness function appears in the exponential. As a result, the description of the proton structure function data would not be as good as with the original fit [34]. Also, as shown in Ref. [32], this parametrization tends to slightly underestimate the coherent γp cross section measured at HERA. In principle one should perform a new fit to the HERA structure function data with the modified IPsat model parametrization. As the purpose of this work is to study the effect of the subnucleon scale fluctuations on ultraperipheral heavy ion collisions this fitting is left for future work. 2 The dipole-nucleus amplitude N A is obtained by using an independent scattering approximation, similar to Refs. [18,33] where b i are the transverse positions of the nucleons in the nucleus. The interpretation here is that 1 − N is the probability not to scatter off an individual nucleon, thus is the probability not to scatter off the entire nucleus.

III. DIFFRACTIVE DEEP INELASTIC SCATTERING
The scattering amplitude for diffractive vector meson production in γ * -nucleus scattering can be written as [38] A γ * A→V A Here the momentum transfer is where P and P are the momenta of the incoming and outgoing nucleus, respectively. The subscripts T and L refer to the transverse and longitudinal polarization of the virtual photon with virtuality Q 2 . In ultraperipheral events, the photons are approximately real (Q 2 = 0) and only the transverse component contributes. The scattering amplitude (5) can be interpreted as follows: first, the incoming virtual photon fluctuates into a quark-antiquark dipole with transverse separation |r|, the quark carrying the momentum fraction z. This splitting is described by the virtual photon wave function Ψ (see e.g. [39]). As discussed previously, the elastic scattering amplitude for the dipole to scatter off the nucleus is N A (r, b, x P ). Finally, the vector meson V is formed, and the qq → V formation is described by the vector meson wave function Ψ V . In this work we use both the Boosted Gaussian and Gaus-LC wave function parametrizations from Ref. [38] in order to estimate the model uncertainty related to the formation of the J/Ψ (see also Ref. [40] for a recent more rigorous calculation of the vector meson wave functions).
Two phenomenological corrections to the diffractive cross sections are included. First, equation (5) is derived assuming that the dipole amplitude is completely real, which makes the diffractive scattering amplitude purely imaginary (in case of a rotationally symmetric dipole amplitude). A correction for the presence of the real part is necessary. Secondly, the skewedness correction that takes into account the fact that in two-gluon exchange processes the gluons in the target are probed at different values of Bjorken-x is also included. These corrections are discussed in more detail in Appendix A.
The coherent diffractive cross section is obtained by averaging the diffractive scattering amplitude over the target configurations and taking the square [13,38]: Here the brackets refer to averages over different configurations of the target. The incoherent cross section is obtained by subtracting the coherent cross section from the total diffractive cross section. It takes the form of a variance of the diffractive scattering amplitude [13] (see also Refs. [14,15,17,18]): The cross sections are related to Fourier transforms of the dipole-nucleus amplitude from coordinate space to momentum space, and the transverse momentum transfer ∆ is the Fourier conjugate to b − (1 − z)r. Here the impact parameter b points to the center of the dipole from the center of the nucleus, and the factor (1 − z)r appears due to the use of non-forward wave functions [38,41]. The dependence on b shows that diffractive vector meson production is sensitive to the impact parameter profile, in contrast to calculations of proton structure functions where the impact parameter integral merely affects the overall normalization. This makes diffractive scattering a sensitive probe of the internal geometric structure of hadrons and nuclei. In particular, larger momentum transfers probe smaller distance scales, which we will show explicitly later. Because the incoherent cross section (7) has the form of a variance, it is sensitive to the amount of fluctuations in coordinate space. In Ref. [18] the incoherent cross section was calculated analytically assuming that the impact parameter dependence of the dipole amplitude factorizes, and is explicitly proportional to e −b 2 /(2Bp) . In that case, the dipole amplitude does not saturate to unity at large dipoles or at large densities. As we calculate the target averages explicitly using a Monte Carlo method (similar to SARTRE [19]), we do not need to rely on this approximation. We note that for the J/Ψ production at the LHC, usage of the factorized approximation for the dipole amplitude reduces the coherent cross section by approximately 15% when no subnucleonic fluctuations are included.
The averages over target configurations are calculated by sampling hundreds configurations. This involves sampling nucleon positions from a Woods-Saxon distribution and the subnucleonic structure for each of the nucleons as described above.
The structure of the target is probed at the scale which can be interpreted as the longitudinal momentum fraction of the nucleon carried by the color-neutral "pomeron". Here W is the center-of-mass energy in the photon-nucleon scattering and m N and M V are the nucleon and the vector meson masses, respectively.
(9) Here y is the rapidity of the vector meson, and the total photon flux n Ai (ω) is obtained by integrating the photon flux of the nucleus over the impact parameter in the region |b| > 2R A . For an explicit expression, the reader is referred to Ref. [42]. The photon flux can be calculated by noticing that the photon energies are ω 1 = (M V /2)e y and ω 2 = (M V /2)e −y . The center of mass energy squared of the photon-nucleon system is W 2 = √ s N N M V e ±y , and the Bjorken-x probed in the Note that there are two different contributions to the vector meson production at y = 0. Either a large-x photon scatters off a small-x gluon, or vice versa. This limits the applicability of the framework at forward and backward rapidities, where a significant contribution to the cross section comes from a process where the center-ofmass energy of the photon-nucleon scattering is small.

V. RESULTS
The coherent and incoherent J/Ψ production cross sections in ultraperipheral Pb+Pb collisions at √ s N N = 5.02 TeV calculated using the Boosted Gaussian wave function are shown in Fig. 1. The cross sections at y = 0 are first calculated without subnucleonic fluctuations, and then including both the geometric and Q s normalization fluctuations for the nucleons. Only the three first coherent peaks are shown, as calculating the average becomes numerically challenging at high |t| [19]. Further, it will be hard to measure the coherent cross section in this regime. We find that the coherent cross section is slightly reduced when the subnucleon scale fluctuations are included. This change is caused by the modification of the impact parameter dependence of the dipole amplitude discussed in Sec. II.
Below |t| 0.25 GeV 2 , where one is sensitive only to fluctuations on length scales larger than the nucleon size, the incoherent cross section is approximately the same with and without nucleon structure fluctuations. On the other hand, at larger |t| subnucleonic fluctuations clearly modify the slope of the incoherent cross section. The |t| value 0.25 GeV 2 , where the incoherent cross section becomes sensitive to subnucleonic fluctuations, corresponds to a distance scale ∼ 0.4 fm, which is of the order of the sizes of the gluonic hot spots in the nucleon (their root mean square radius is 2B q ≈ 0.24 fm). Preliminary ALICE data on exclusive J/Ψ production [43] show that the change in slope occurs around √ −t ≈ p T ∼ 0.5 GeV, which is in quantitative agreement with our results. Next, we compare our results with the total (t integrated) J/Ψ production cross sections measured by AL-ICE [8,9] and CMS [10] at √ s N N = 2.76 TeV as a function of the J/Ψ rapidity. The results for coherent and incoherent diffraction are shown in Figs. 2 and 3, respectively. We show results obtained by using both the Boosted Gaussian and Gaus-LC wave functions for the J/Ψ. Similarly to previous literature [27], we find that the different wave functions mainly affect the overall normalization of the cross section. In Fig. 2 we see that the coherent cross section is somewhat reduced when subnucleonic fluctuations are included. This change is comparable to the model uncertainties related to the J/Ψ wave function. In particular, we find that replacing the Boosted Gaussian by the Gaus-LC parametrization, both the coherent and incoherent (shown in Fig. 3) cross sections are reduced by approximately 20%.
Using the Boosted Gaussian wave function, without subnucleon scale fluctuations, the coherent cross section is significantly overestimated, as in Ref. [27], and the incoherent cross section agrees with the data. Including geometric fluctuations increases the incoherent cross section almost by a factor of 2, and the results are consistently above the data for both processes. The rapidity dependence of the data is well reproduced. Results with subnucleon scale fluctuations obtained with the Gaus-LC wave function are close to the experimental data for the coherent cross section and slightly higher for the incoherent cross section. At midrapidity, the ALICE coherent cross section datapoint is overestimated by 1.3σ and the incoherent cross section by 2σ. We note that there is tension between the HERA e + p data where the coherent cross section is underestimated by our model with subnucleonic fluctuations [32].
To reduce model uncertainties related to the vector meson wave function we study the ratio of incoherent to coherent cross section. Our results are compared to the ALICE data [8] at √ s N N = 2.76 TeV in Fig. 4. It can be seen that inclusion of the subnucleonic fluctuations brings this ratio to a level compatible with the experimental data. The ratio is found to be approximately independent of rapidity and we confirm that the dependence on the vector meson wave function is very weak. Predictions for the J/Ψ rapidity distribution in ultraperipheral Pb+Pb collisions at √ s N N = 5.02 TeV are shown in Figs. 5 and 6. As for the lower energy, we find that the subnucleon scale geometric fluctuations have a large effect on the incoherent cross section. The ratio of the cross sections for the two processes for √ s N N = 5.02 TeV decreases due to larger saturation effects on the incoherent cross section [18]. Numerically this effect is small in the LHC energy range, and we find that at √ s N N = 5.02 TeV the ratio decreases by 0.5 . . . 3%.
Finally, in Fig. 7 we show |t|-spectra for J/Ψ pro- GeV Au+Au collisions at RHIC, corresponding to W = 25 GeV. This corresponds to x P = 0.015, which is at the edge of the validity of our model. Especially the skewedness and real part corrections together are almost 100%, which makes the absolute normalization unreliable (see Appendix A and Fig. 8). The spectra are calculated using the Boosted Gaussian wave function. Integrating the cross sections over t, we get 106 µb for the coherent and 62 µb for the incoherent cross section with subnucleon fluctuations (with Gaus-LC wave function the cross sections are 100 µb and 52 µb). The corresponding cross The Bjorken-x evolution in this work comes directly from the Q s evolution in the IPsat model. Thus, the amount of fluctuations and the size of the hot spots do not change as a function of x or center-of-mass energy. In principle the characteristic length scales (∼ Q −1 s (x)) depend on the energy and recent explicit calculations show that protons grow and fluctuations are reduced when Bjorken-x is decreased [21,45]. If that is the case, we would expect the incoherent cross section to grow more slowly with energy than in our calculation.

VI. CONCLUSIONS
We have shown that the subnucleon scale fluctuations, in particular geometric fluctuations of the nucleon shape, contribute significantly to the incoherent J/Ψ production cross section at |t| 0.25 GeV 2 , measured in ultraperipheral heavy ion collisions at the LHC and RHIC. This is the first main result of this work: we expect to see the incoherent t slope of the J/Ψ production cross section to change at the value of t corresponding to the distance scale of the subnucleonic fluctuations. Compared to the case where the only contribution to the fluctuations originates from fluctuating nucleon positions, the |t| integrated incoherent cross section increases almost by a factor of 2 when geometric fluctuations of the nucleon shape are included.
When comparing to experimental data, the cleanest presented result is the ratio of the incoherent and coherent cross sections, which eliminates a large part of the model uncertainties. It increases by a factor of two when subnucleonic fluctuations are included. This is the second main result of this work, and explains previous tension between the employed dipole model and experimental data in Ref. [27], where one generally overestimated the coherent and underestimated the incoherent cross section.
In this work the energy evolution only affects the saturation scale of the nucleons. Our framework does not involve any possible smoothening of the nucleon or nucleus as one evolves towards small-x, discussed e.g. in Refs. [21,45]. As shown in [21], one could expect that the smoothening of the proton slows down the growth of the incoherent cross section with energy. Including these effects in our calculation explicitly by solving JIMWLK evolution equations [46][47][48][49] as done in e.g. [50], is left for future work. effect in photon-nucleus scattering. This correction is taken into account by multiplying the obtained cross sections by a factor 1 + β 2 .
The second phenomenological correction we include the skewedness correction, which takes into account the fact that in the two-gluon exchange the gluons in the target are probed at different longitudinal momentum fractions x 1 x 2 ≈ x P [51][52][53]. In the IPsat model the collinear factorization gluon distribution x P g(x P , µ 2 ) is corrected to correspond to the off-diagonal (or skewed) distribution, which depends on both x 1 and x 2 , by multiplying it by a skewedness factor R g . Following the prescription of Ref. [38] we get R g = 2 2λg+3 Γ(λ g + 5/2) √ π Γ(λ g + 4) (A3) with λ g = d ln x P g(x P , µ 2 ) d ln 1/x P .
For photon-nucleus scattering the skewedness correction is approximated by calculating its effect on the photonproton scattering, and using the obtained correction factor.
Especially the skewedness correction is numerically important and needed to describe the HERA diffractive measurements. The effect of real part and skewedness corrections on the total coherent diffractive cross section is shown in Fig. 8. When corrections are calculated, no proton fluctuations are taken into account. The correction grows rapidly towards small rapidities (small y). Note that when J/Ψ production is calculated at nonzero rapidity, there are always large-x contributions (corresponding to negative y) and small-x contributions (positive y) to the cross section, see Eq. (9).