Investigating saturation effects in ultraperipheral collisions at the LHC with the color dipole model

We investigate saturation effects in ep scattering as well as in ultraperipheral pA and AA collisions at small x with four variants of the impact parameter dependent color dipole model: with and without gluon saturation and with and without a novel mechanism that suppresses unphysical dipole radii above the confinement scale, a problem not addressed by most implementations. We show that ep scattering at HERA can be very well described by any of the four variants. When going from ep to eA scattering, saturation effects are expected to increase as ∼A1/3. In lieu of an electron-ion collider, we confront the different versions of the dipole model with data recorded in ultraperipheral collisions at the LHC in order to estimate the sensitivity of the data to gluon saturation in the target nuclei. We find that ultraperipheral PbPb collisions indicate strong saturation effects while pPb collisions turn out to not have any discriminating power to distinguish saturation from non-saturation scenarios.


Introduction
Over the last two decades the Color Dipole Picture [1][2][3][4][5][6] has been developed to study high energy scattering in QCD. An attractive feature of the dipole approach to high-energy interactions is that it gives a clear interpretation of the physics at small values of Bjorken x. Even in its simplest form, it turns out to be rather successful in describing the total [1,7,8] and diffractive [2,9,10] virtual photon-nucleon cross-sections. Its generalizations are now commonly used for the parametrization of data from HERA [6,11,12]. In the scattering of leptons off hadrons at high energy, high-density phases of partons are created and non-linear effects become important. In this "saturation" regime, the basic properties of perturbative QCD, such as collinear factorization and linear evolution, break down. This unique form of strongly interacting matter is called a Color Glass Condensate (CGC) [13]. The color dipole model has been originally created to study scattering in this regime.
From a phenomenological point of view, the key is to connect the experimental probes to the scattering of a dipole. In the case of lepton-hadron collisions this is obvious: a lepton undergoes hadronic interactions via a virtual photon and one can regard the interaction as the fluctuation of the virtual photon into a quark-antiquark pair which then interacts. In hadron-hadron collisions however, there are no such virtual-photon probes. The exception are ultraperipheral collisions (UPC) where, at very large impact-parameters between the colliding hadrons, the long range electromagnetic force becomes dominant over short-range QCD. The intensity of the electromagnetic field, and therefore the number of photons in the cloud surrounding the nucleus, is proportional to the square of the hadron's electric charge. Thus these types of interactions are highly favored when heavy ions collide. UPC are currently extensively measured at RHIC [14][15][16][17] and LHC [18][19][20][21][22][23][24][25][26].
Various attempts have been made to describe UPC data in the dipole model to test if saturation effects are present [27][28][29][30][31][32][33][34][35]. The parameters that provide the necessary non-perturbative input to these models are commonly determined through fits to high-precision structure functions measured at HERA in ep collisions [4][5][6]. So far, the comparison of color dipole models with UPC data has provided no clear evidence for or against saturation. Failure (or success) in describing UPC data with the dipole model does not imply the absence (or presence) of saturation but could also point to shortcoming of the model itself. Improvements, such as NLO calculation, are under way [36,37]. Another issue, one we address in this paper, is that the dipole model allows for large dipole radii beyond the confinement scale which are unphysical. We will show how this affects the comparison with data using the diffractive event generator Sartre [38,39]. To overcome the issue of unphysical large radii we implemented a damping mechanism inspired from what was originally introduced in [40]. To do so, new dipole parameters had to be obtained by fitting the modified dipole model to HERA data.
In this paper we also consider two versions of the dipole model, one with saturation (bSat 1 ) and one without (bNonSat). The latter is obtained by linearizing the dipole cross-section. This comparison allows us to test the actual sensitivity of UPC data to saturation effects.

Nonlinear effects in the Dipole Model
In Deep Inelastic Scattering (DIS) electron-hadron scattering a virtual photon interacts electromagnetically with a parton in the hadron. At small parton momentum fractions x, this process can be seen as the virtual photon fluctuating into a quark anti-quark color dipole which subsequently interacts via one or more gluon exchanges with the proton. This process is described by impact parameter dependent dipole model [3,4]. The DIS process has the following crosssection: This can be seen as a three part process where the virtual photon splits into a color dipole, which interacts with the hadron and then recombines. The splitting and recombination of the virtual 1 In literature, bSat (bNonSat) is also known under the name IPSat (IPNonSat).
2 photon is described by the transversely and longitudinally polarized wave functions: Here, z is the quark's momentum fraction of the photon, m f the quark mass, Q 2 the photon virtuality, x the gluon's momentum fraction of the hadron, and g is the DGLAP longitudinal gluon density [41][42][43][44]. The interaction between the dipole and the hadron is described by the dipole cross-section, which can be seen as the exchange of one gluon [3]: There are three types of phenomenological modifications in the dipole cross-section which address non-linear effects. Firstly, the proton thickness T p (b) cannot be calculated from first principle. However, for exclusive diffraction, the impact parameter b is an observable, and it can be experimentally accessed via its Fourier conjugate, the Mandelstam variable t. Investigations of the t distributions in exclusive J/ψ production have shown that the thickness is well described by a Gaussian T p (b) = exp(−b 2 /2B p )/2πB p , with the parameter B p = 4 GeV −2 [4][5][6] . The proton's thickness will not be the focus of this paper. It is further expected that the dipole size r is dampened for radii larger than the confinement scale. There is an inherent damping effect in the Bessel functions in the wave overlap, since at large r these are suppressed as ∼ exp(−r ), giving an effective size of the dipole of r ∼ 1/ . If is small, which can happen at small Q 2 and at z 0, 1, this does not give enough suppression of large dipoles. The authors in [3,4] solve this by introducing artificially large masses of the light quarks similar to pion or ρ-meson masses ensuring a sizable analogous with a meson enhancement. Later approaches [5,6] have ignored this problem and allowed for any , and thereby allowed for unnaturally large dipoles. A novel approach is to add a Gaussian suppression to the dipole cross-section by hand, inspired by Flensburg et al. [40]: where r soft is the modified dipole size used in the dipole cross-section, r pert is the dipole size given by the wave overlap. For small dipoles, r soft ≈ r pert as expected. Here R shrink is a free parameter. This is the approach we will adopt in this paper. It is worth noting that there are two effects of damping that have the opposite effect on the dipole cross-section. The shift r pert → r soft , where r soft < r pert , causes a trivial, x-independent suppression due to the r 2 factor in dσ dip /d 2 b. However, the shift in r also affects xg(x, µ 2 ) since µ 2 = C/r 2 + µ 2 0 . The increase in the scale, ∼ 1/r 2 soft , increases typically the gluon density thus increasing the cross-section. This introduces also an x dependence to the damping.
The third consideration is related to saturation, which also suppresses large dipoles as well as very small values of x. Saturation effects are important in the region of phase-space where the gluon wave-functions are overlapping, and large dipoles have an increased amplitude for 3 interacting with multiple gluons and therefore pick up their correlations. Also, at small x, the gluon wave-functions becomes spatially larger and therefore overlap more. The saturated version of the dipole model may in principle be derived from the Color Glass Condensate effective theory for QCD [45], and is of the following form: At small r, this expression becomes equal to Eq. (2). At first glance, this too seems to give a Gaussian suppression for large dipoles. However, the r-dependence in the scale of the DGLAP gluon density makes the suppression more involved. In Eq. (4) large gluon densities are also suppressed as expected. Eq.(4) is referred to as the bSat model, while Eq. (2) is referred to as the bNonSat model. In order to disentangle large dipole effects we believe that we need to have high precision measurements of the gluon density in the nucleus. When going from ep to eA collisions, the confinement scale is expected to remain unchanged while the saturation scale gets an "oomph" of ∼ A 1/3 . Thus, by investigating the A dependence of the cross-section one should be able to disentangle these non-linear effects. This will be one of the main foci of the future Electron-Ion Collider (EIC) [46] which will be constructed in the US in the next decade. In lieu of an EIC, ultraperipheral collisions measured at RHIC and LHC offer the possibility to study similar processes in photoproduction, where Q 2 ≈ 0. In this paper we will attempt to disentangle the large dipole effects using combined e ± p DIS data from HERA I+II as well as pPb and PbPb UPC data for exclusive vector meson production at the LHC. For this purpose we have implemented the calculations in the event generator Sartre [38,39].

Fits to inclusive DIS at HERA
The proton structure functions F 2 (x, Q 2 ) and F L (x, Q 2 ) can be written in terms of the total photon-proton cross-sections as: Following [6], we take the DGLAP gluon density at a starting scale µ 2 0 = 1.1 GeV 2 to be xg(x, 6 , using variable flavour scheme when evaluating the strong coupling α s and solving the DGLAP evolution for the gluon density. The strong coupling satisfies For heavy quarks, the Bjorken x is replaced by We fit to the combined reduced cross section data from HERA I and HERA II [47,48] and the combined reduced cross section for charm data [49,50] from the H1 and ZEUS experiments. The reduced cross section is given by: We include data in the range x < 0.01, or x f < 0.01, and 1.5 ≤ Q 2 ≤ 50 GeV 2 . This gives us 409 data points for σ r and 34 points for σ cc r . The results from the fits including confinement are 4 Model  [47][48][49][50]. Here B p = 4 GeV −2 , µ 2 0 = 1.1 GeV 2 , m b = 4.75 GeV, and m t = 175 GeV. X+Y points means X points for inclusive and Y points for charmed reduced cross section. The table also contains the results from the previous fit [6] for reference. For inclusive DIS, bNonSat with damping has χ 2 /Ndf=1.073, and bSat with damping has χ 2 /Ndf=1.262.
presented in Table 1 as well as the result of the previous fit [6] for comparison. When selecting the data points included in the fit, we used m c = 1.321 GeV, which differs slightly from the choice in [6]. This causes a negligible difference in the number of data points used for our fit.
We observe that the bNonSat model with damping fits the HERA data slightly better than the bSat model, although all models describe the data reasonably well. We also see that adding damping effects to the bSat model does not significantly alter the resulting fit quality. For fits with the explicit damping model, the light quark masses are substantially smaller, as expected. Also, for bNonSat, adding damping allows for a slower growth of the gluon density for small x, as seen by the values of λ g . Our fit results indicate that there are no definitive hints for saturation effects in the HERA data, since the data can be equally well described using a Gaussian suppression of large dipoles. Since the fit procedure is identical to that used in [6], we refer to this reference for an in-depth discussion of the sensitivity of the result on assumed quantities such as proton profile width B p , bottom quark mass m b , starting scale for DGLAP evolution µ 0 , and the choice of wave function.
The resulting dipole cross-sections are shown in Fig. 1 (a) and (b) for proton and lead respectively, illustrating how the damping suppress large dipoles in the bNonSat model. We also note that adding damping to the bSat model increases the dipole cross-section slightly at smaller r. To anticipate the effects of damping and saturation in exclusive J/ψ production in UPC, we also show in (c) the dipole cross-sections multiplied by the wave-function overlap between the incoming virtual photon and the produced J/ψ meson, which by itself is also suppressing dipole radii larger than the typical size of the vector meson r ∼ 0.06 fm. We are using the "Boosted Gaussian" wave overlap parametrization from [4]. We note that for lead, both saturation and damping have significant impact, while for the proton the damping mechanism completely takes away any saturation effects. In Fig. 1 (d) we show how damping is affecting the saturated dipole cross-section for protons at different impact parameters, and we see that the damping effect becomes stronger at larger impact parameters.

Exclusive Vector Mesons in UPC with Sartre
In ultraperipheral collisions, the interacting hadrons are so far apart that the long range electromagnetic force dominates over the short range strong force. These interactions are therefore very similar to DIS with exchange of pseudo-real photons with Q 2 ≈ 0. For exclusive diffractive production of vector mesons in photoproduction the differential cross-section can be written as:  6 where y = ln(2E γ / m 2 V + p 2 V⊥ ), giving dE γ /dy = E γ and x I P is the momentum fraction of the exchanged gluons with respect to the probed hadron. Here n γ is the flux of photons interacting with the proton. In the derivation of the amplitude only the imaginary part is taken into account. The correction for the missing real part is given by β = tan(λ · π/2), where λ = −∂ ln A T,L /∂ ln x I P [4] . The above cross-section treats the diffractive exchange as one or multiple two-gluon exchanges where the two gluons shield each other's colors, ensuring that the interaction does not change the proton's quantum numbers. However, to account for the possibility of the two gluons having different momentum fractions x I P , the amplitude is multiplied by a skewedness correction R g given by: where λ skew = − ∂ ln(xg(x,µ 2 )) ∂ ln x . As the skewedness corrections are formally needed only in the case of linear DGLAP gluon densities, and the dipole model modifies this density significantly, it is not theoretically clear whether or not it should be applied in this case. However, it has historically been required to describe J/ψ production in ep collision at HERA. In the following we will apply the skewedness correction to the proton dipole, where it is better motivated, but omit it for the ion dipoles. Both corrections become smaller with decreasing x. In the range relevant for this work, the real part correction is in the range 1 + β 2 = [1.17, 1.53] while for the skewedness correction R g = [1.5, 2.0]. In the case of γ * A scattering, the bSat and bNonSat dipole cross-sections become: where b i are the positions of the nucleons in transverse space for a given nucleon configuration. The total cross-section is then given by: where the average is taken over initial state nucleon configurations. The coherent part of the cross-section, in which the struck nucleus stays intact after the interaction, is given by: The incoherent case, where the struck nucleus becomes excited in the interaction and subsequently de-excites by emitting a photon, one or many nucleons, or fragments, is given by the difference between the total and coherent cross-sections. In Sartre we generate 500 nuclear configurations and average the amplitude over these. This has been shown to be sufficient for the cross-sections to converge [39]. For the photon flux we follow the model used in the STARLIGHT generator [51]. We have implemented these processes in the Sartre event generator which give an exclusive final state with exact four-momenta of all incoming, intermediate and final state particles. With Damping Figure 2: A comparison between Sartre and measured exclusive J/ψ data from ALICE [20,21] in pPb UPC at √ s = 5.02 TeV. We show the results using dipole models without (with) damping of large dipole radii on the left (right) hand side.
In Fig. 2 we show cross-sections from two measurements in pPb collisions by the ALICE collaboration [20,21] with √ s = 5.02 TeV and compare them with results from the Sartre event generator. These measurements cover a rapidity range of |y| < 4 which corresponds to 10 −5 ≤ x I P ≤ 10 −2 . In the left panel we show the Sartre results without damping, and in the right panel with damping, both in the same kinematic range as the data and integrated over t. We see that both versions of bSat are able to describe the data well. The bNonSat model without damping (left panel) deviates slightly more from bSat at all rapidities in the range, while the damped version of bNonSat (right panel) does not differ significantly from bSat. However, any differences are within the uncertainty range of the data.
In Fig. 3 we show the results from events generated for ultraperipheral PbPb collisions at √ s = 2.76 TeV, and compare them with measurements from ALICE [52] and CMS [24]. These measurements cover a rapidity range of |y| < 3.5 which corresponds to 3 · 10 −5 ≤ x I P ≤ 3 · 10 −2 . Data and Sartre are integrated over t. For symmetrical beams, it is not experimentally possible to distinguish if the photon comes from beam 1 or beam 2, and the total cross-section is: Therefore, away from y = 0 (x I P = 10 −3 ) the cross-section is a linear combination of larger and smaller x I P . The dipole model is only applicable for x I P 0.01, which we relax a little since in this case it also imposes a lower bound on x I P . Sartre is therefore only able to access |y| < 2.9 at this beam energy configuration. Both bSat versions can describe the coherent data well. Sartre appears to favor a non-saturated model for the incoherent data as the bSat description lies below the data. However, this is to be expected, since the incoherent cross-section is proportional to gluon fluctuations, and Sartre only includes fluctuations of the initial nucleon configuration in its model. We expect that this description will change once we have also included subnucleonic fluc-8  [52] and CMS [24] of PbPb UPC at the LHC at √ s = 2.76 TeV. The left (right) side show the results using dipole models without (with) damping. The arrows on the incoherent cross sections indicate where the curves would be if we included subnucleonic gluon fluctuations in the calculations, which according to [31] contributes a factor 1. 8. tuations and saturation scale fluctuations in our model following the prescription by Mäntysaari and Schenke [31,53,54]. In [31] the authors show that subnucleonic fluctuations increases the incoherent cross section in this beam configuration by a factor of 1.8. This is indicated by arrows in Fig. 3. We see that even with damping there is a large difference between bSat and bNonSat, at y = 0 which corresponds to x I P = 10 −3 , as anticipated from Fig. 1 (c). This difference is even larger in the incoherent cross-section.
In Fig. 4 (a) we show the non-saturated to saturated cross-section ratios as a function of x I P in exclusive J/ψ production for incoherent ePb and coherent ep and ePb collisions. We see that the largest saturation effects occur in the incoherent ePb cross-sections, while for ep there is little sign of saturation, especially in damping mode. In Fig. 4 (b) we depict the same ratio as a function of t for 10 −5 ≤ x I P ≤ 10 −2 . Note that the coherent cross-section in the bSat models decreases faster as a function of |t| than in the bNonSat models. We further note that at larger |t|, where the incoherent part of the cross section dominates, the ratio becomes large. For the proton, the ratio remains near unity for the entire t-spectrum. The large difference between bNonSat and bSat in Fig. 3 comes from the integral over all t. Fig. 4 (c) illustrates the ratio of cross-section without dipole-radius damping over that with damping as a function of x I P . We see that the x I P dependence of the damping is small and verify that it is independent of the nuclear species. The damping effects are identical for coherent and incoherent scattering. Naively the damping is not expected to have any dependence on x I P . However, the absence of damping in the dipole cross section is compensated for by a smaller value of λ g in the fit, causing a slower growth in the gluon density at small x I P , which can be seen in Table 1. There is also a direct dependence on the dipole radius, and therefore on the damping, in the factorization and renormalization scales µ 2 = C/r 2 + µ 2 0 , which gives a larger effect in the DGLAP evolved gluon density at small x I P and large µ 2 . This is further enhanced by the J/ψ  wave function overlap. What is seen in Fig. 4 is the combination of these effects. In Fig. 5 we show a comparison of the W-dependence of the total cross-section for exclusive J/ψ photoproduction, γ * + p → J/ψ + p, between our four model variants and measurements from the H1 [55,56], ZEUS [57] and ALICE [20,21] collaborations. In addition we depict our predictions for γ * + p → J/ψ + Pb, which will be tested at a future EIC. Again, we see that the difference between the models in γ * p interactions is not large enough for the data to distinguish between the saturation and non-saturation scenario, while in γ * Pb there is a pronounced separation between the two.

Conclusions
We have introduced a damping for dipole radii larger than the confinement scale in the bSat and bNonSat models to improve the color dipole model. We show that this new model with newly derived parameters can describe HERA data well. The dipole model is implemented in the Sartre event generator. Our studies demonstrate that the damping improves the description of the non-saturated dipole model for ultraperipheral pPb collisions at the LHC. However, the precision of the current data does not allow to discriminate between saturation and non-saturation scenarios for pPb collisions.
The improved model describes ultraperipheral PbPb collisions quite well. Saturation effects in coherent PbPb collisions appear to be significant for all rapidities. A comparison of data from incoherent interactions will need to be further improved to include gluon fluctuations in the nucleons.  Figure 5: A comparison of the W dependence of the different models for γ * + p → J/ψ + p and γ * + Pb → J/ψ + Pb. The proton case is compared to measurements from H1 [55,56], ZEUS [57] and ALICE [20,21] collaborations. All cross-sections are integrated over t.