Systematics of strong nuclear amplification of gluon saturation from exclusive vector meson production in high energy electron-nucleus collisions

We show that gluon saturation gives rise to a strong modification of the scaling in both the nuclear mass number $A$ and the virtuality $Q^2$ of the vector meson production cross-section in exclusive deep-inelastic scattering off nuclei. We present qualitative analytic expressions for how the scaling exponents are modified as well as quantitative predictions that can be tested at an Electron-Ion Collider.


I. INTRODUCTION
A striking observation made by the HERA deep inelastic scattering (DIS) experiments on electron-proton collisions is the rapid growth of the gluon densities at high energies or, equivalently, at small values of longitudinal momentum fraction x carried by the gluons [1,2]. Weak coupling QCD studies predict that when gluon occupancies become of order 1/α s , where α s is the QCD coupling constant, the perturbative bremsstrahlung of gluons is balanced by strong repulsive gluon self-interactions. This dynamics generates a phenomenon called gluon saturation that tempers the growth in gluon distributions [3,4]. For every Q 2 Λ QCD , with Λ QCD representing the intrinsic QCD scale, there is an x value for which the maximal occupancy is reached. Equivalently, at small x, for each value of x there exists saturation scale Q 2 s,p (x) in the proton; transverse momentum modes below this scale have maximal occupancy. High momentum modes k ⊥ Q s,p (x) asymptote to the usual perturbative QCD dynamics. The underlying physics of the onset of and the many-body dynamics of gluon saturation is captured in a Color Glass Condensate (CGC) effective theory [5][6][7][8][9][10] which allows for efficient computations.
Despite the large amount of precise HERA data, the evidence for saturation effects is unclear even though recent analyses increasingly point in this direction [11,12]. One reason for this uncertainty is that the proton's saturation scale Q 2 s,p (x) is not very large at the x values probed in the HERA experiments. Discussions of gluon saturation in hadron-hadron collisions can be found in [13,14]. The search for gluon saturation in particular and the many-body non-linear dynamics of gluons in general is a major motivation for the Electron Ion Collider proposal [15] in the US and the LHeC proposal at CERN [16].
The physics of gluon saturation is universal. All nonperturbative information specific to the quantum numbers of a given nucleus is encoded in the initial conditions; at asymptotically small x, at a fixed impact parameter, an external probe would be unable to distinguish between a heavy nuclear target and a proton. All memory of their initial conditions would be lost and their saturation scales would be nearly identical [17]. However for realistic values of x, as shown explicitly in the McLerran-Venugopalan (MV) model [5][6][7], nuclear information is contained in the atomic number (A)dependence of the saturation scale. While not universal in the larger sense of complete independence of the dynamics from the initial conditions with which the system is "prepared", it is nonetheless remarkable that the dynamics of a nucleus at high energies is controlled primarily by a single A-dependent scale. In the MV model, simple kinematic and dynamical arguments suggest that ; this expection is borne out in detailed models [18] that we shall discuss further shortly.
The search for gluon saturation will require that we identify measurements where its onset will show strikingly different systematics from those seen in its absence. Based on the arguments outlined, it is widely believed that effects of gluon saturation are likely precocious in high-energy (e + A) DIS off nuclei. We will demonstrate here that the onset of gluon saturation is especially dramatic in exclusive vector production off large nuclei. The simple reason why exclusive processes are attractive measurements to pursue is that the perturbative QCD crosssection in such processes is proportional to the nuclear gluon distribution squared [19]. Since the gluon distribution grows rapidly at small x, such measurements are especially sensitive when gluon distributions saturate.
We will discuss here the Q 2 and A dependence of exclusive vector meson production in e+A collisions within the CGC framework. Exclusive production of vector mesons in this saturation picture of e + A DIS has been studied previously in the literature [20][21][22][23][24]. A recent development in exclusive vector meson electroproduction is the extraction of event-by-event fluctuations of gluon spatial distributions in the proton from HERA data [25]. These results are very relevant to understanding recently discovered "ridge-like" multiparticle correlations in protonproton and proton-nucleus collisions [26,27] and can also be studied in future in e + A collisions [24,28]. There has been recent progress in developing the theory beyond leading log x accuracy [29]; we expect that these developments will in future help us refine our estimates.
Our focus here will be to demonstrate that gluon sat-uration gives rise to very specific (and strong) dependencies in the A and Q 2 scaling of exclusive vector meson cross-sections for Q 2 > Q 2 s,A that are qualitatively different from those for Q 2 < Q 2 s,A , where Q s stands for the saturation scale of the target. Exclusive processes off nuclei are currently being studied in in ultraperipheral heavy-ion collisions at RHIC and the LHC where one nuclei acts as a source of (quasi) real photons [30,31] for the other; there are several recent works that explore the role of small-x dynamics in these measurements [32][33][34][35][36][37][38]. However in these reactions, the photon virtuality is limited to be Q 2 ∼ 1/R 2 A , where R A is the nuclear radius. Further, studying the A dependence by varying the heavy-ion beams is experimentally challenging.
These limitations do not exist at the EIC [15] where exclusive processes can be measured over a wide range of Q 2 , A as well as rapidity separation y P = ln(1/x P ) between the vector meson and the target. The variable x P is the equivalent of Björken x for an exclusive process; it has the parton model interpretation of being the momentum fraction of a parton within the colorless ("Pomeron") exchanged between the virtual photon probe and the nuclear target. The kinematical coverage of the EIC is discussed in detail in Ref. [39].

II. EXCLUSIVE VECTOR MESON PRODUCTION
At high energies, DIS can be formulated conveniently in a dipole picture, whereby the incoming photon with virtuality Q 2 splits into a quark-antiquark dipole and the dipole subsequently scatters elastically off the target with the amplitude N ; in exclusive vector meson production, the dipole hadronizes into the vector meson. The virtual photon to quark-antiquark dipole splitting amplitude Ψ γ * →qq can be computed in quantum electrodynamics [9]. The dipole to vector meson hadronization amplitude Ψ qq→V M is non-perturbative. We therefore have to rely on phenomenological parametrizations of the vector meson wavefunction; we will employ here the Gauss-Light Cone parametrization from [21]. Because the time scales for hadronization are far greater than the time scale of the interaction of the dipole with the target, ratios of the exclusive vector meson production crosssections in different nuclei should at most weakly depend on the the details of these non-perturbative wavefunctions.
The scattering amplitude for exclusive vector meson production can be written as [21] Here z is the fraction of the momentum of the photon carried by the quark (the antiquark carries the remain-ing fraction 1 − z), ∆ ≈ √ −t is the transverse momentum transfer from the target and σ dip is the dipoletarget cross-section. The difference between forward and non-forward wavefunctions gives the extra factor [21,40]. We anticipate that this factor has only a small effect at small x and its role will be further diminished in the ratios we will examine and will not include it in our numerical calculations. The coherent diffractive vector meson production cross-section reads as where V denotes the vector meson of interest. In our numerical computations in the following sections, we will use the IPsat model to describe the dipoleproton cross-section σ dip . The IPsat model contains an eikonalized DGLAP-evolved gluon distribution function (see also Ref. [41]), which guarantees the correct perturbative limit in the dilute region (small x and large Q 2 ) and preserves unitarity because dσ dip /d 2 b T → 2 for large dipoles. The impact parameter dependence is included by having the saturation scale depend on the transverse distance from the center of the proton. In the IPsat model, the expression for the dipole cross-section is [42] dσ dip with the Gaussian proton transverse profile function: The free parameters in the model are the proton size B p measured by the dipole probe and the initial conditions for the DGLAP evolution of the gluon distribution function x g(x, Q 2 ). These were fixed in Ref. [11] by fitting the HERA inclusive DIS data. The dipole amplitude obtained was applied successfully previously in several phenomenological works [36,43,44]. The dipoleproton cross-section can be generalized to nuclei as in Refs. [18,23,42]. It is employed for instance to construct the IP-Glasma initial conditions for A + A collisions [45][46][47]. We also include the so-called skewedness and real part corrections as in Ref. [23]. These corrections are phenemenologically necessary in order to describe the HERA data 1 . These corrections do not affect the A scaling in our discussion. They have a small effect on the Q 2 scaling which mostly cancels in the cross-section ratios.
The generalization of the IPsat dipole cross-section to nuclear targets, following the procedure suggested in [18,42], is given by where σ dip is the total dipole-proton cross-section integrated over the impact parameter. The nuclear transverse density profile T A is obtained by integrating the Woods-Saxon distribution over the longitudinal coordinate, and is normalized as We will supplement our numerical study by simple analytical estimates that capture the underlying physics behind the Q 2 and A scaling of the vector meson crosssections. These analytical estimates are easier to realize in a simpler dipole model that is sufficient for this purpose, the Golec-Biernat-Wusthoff (GBW) model [49]. In this model, the dipole-proton scattering amplitude can be written as and the saturation scale Q s,p implicitly depends on x P . The dipole-nucleus cross-section is obtained from Eq.
In the following, we will present numerical results computed in the IPsat model for J/Ψ and ρ production in the kinematics relevant for the EIC and LHeC as well as analytical estimates in asymptotic kinematics computed in the GBW model. An advantage of studying exclusive ρ electroproduction is that, at moderate Q 2 , the contribution from the saturation region is enhanced relative to the case of the heavier J/Ψ meson. On the other hand, since the ρ is much lighter than the J/Ψ, the lack of a large momentum scale makes perturbative computations questionable at smaller Q 2 values.
We will consider exclusive vector meson production in the high Q 2 region alone in this section and will take up the case of low Q 2 in the following section. In this high virtuality region, Q 2 is the only relevant scale and the mass difference between the J/Ψ and the ρ is less relevant.

A. A scaling
It is instructive to first consider the forward limit of t = 0. Here, in the dilute ("pQCD") region of Q 2 s,A Q 2 , we can expand all exponents in and (6) and get where Q 2 s,A ∼ A 1/3 Q 2 s,p . Performing the b T integral, one gets another factor of A 2/3 from the area, which gives Squaring the amplitude to obtain coherent cross-section, one obtains This scaling for the IPsat model is demonstrated numerically in Fig. 1 which shows the normalized ratio of the exclusive ρ meson cross-sections at t = 0 for Gold and Iron nuclei respectively. Such ratios are desirable because as noted previously model dependencies cancel; this is also the case for systematic uncertainties in the experiments. Fig. 1 demonstrates that for longitudinal polarization, the asymptotic A 2 scaling is seen for Q 2 = 10 3 GeV 2 , though it is within 5% already for Q 2 = 10 2 GeV 2 . At smaller Q 2 , significant suppression seen is due to gluon saturation. The suppression is greater for transversely polarized photons which are more sensitive to larger dipole sizes. The energy or x P dependence of the ratio is weak, with the suppression increasing only slightly when x P is decreased from 0.01 to 0.001. For comparison, we note the values for the saturation scales 2 of Gold and Iron nuclei at b = 0 extracted from the IPsat model generalized to nuclei as shown in Eq. (5). For Gold, we get the value of the quark saturation scale to be Q Since the width of the coherent peak is t max ∼ 1/R 2 A ∼ A −2/3 , the total exclusive vector meson cross-section (integrated over t) scales like The numerical result for this scaling in the IPsat model is shown in Fig. 2. At high Q 2 , the normalized ratio does not go to unity. This is a consequence of oversimplifying the t integral to include just the width of the coherent peak to give the factor A −2/3 , and is valid only at asymtotically large A.
For a given choice of the vector meson wavefunction (longitudinal Gauss-LC from Ref. [21]), the overlap of 2 The quark saturation scale Q 2 s is defined as  this wavefunction with the longitudinally polarized photon wavefunction has the dependence with ε = Q 2 z(1 − z) + m 2 q ≈ Q for Q 2 Q 2 s and |r T | = r. The scalar part of the vector meson wavefunction φ L ∼ z(1 − z)e −r 2 M 2 V limits contributions from dipoles larger than the inverse vector meson mass M V . A stronger limit on dipole sizes, r 1/Q, is set by the Bessel function, and the Q 2 scaling becomes Thus the exclusive longitudinal ρ cross-section (both the total and its value at t = 0) has the Q 2 scaling [50] The numerical result for the dependence of the ρ and J/Ψ exclusive cross-section on Q 2 (scaled by Q 6 ) is shown in Fig. 3. The Q −6 behaviour of the cross-section at high Q 2 is apparent since the scaled cross-section is flat in the region Q 2 10 2 GeV 2 . The A 2 scaling is also visible, as the curves corresponding to different nuclei lie on top of each other in the Q 2 range plotted in the figure.
For transversely polarized photons, we can perform a similar analysis as shown for the longitudinal polarization case in Eq. (11). Neglecting terms proportional to the light quark mass, the wavefunction overlap becomes where the scalar part of the vector meson wavefunction now behaves as φ T (r, z) ∼ z 2 (1 − z) 2 e −r 2 M 2 V . Thus the difference to the longitudinal polarization case corresponds to an extra power of r in the scattering amplitude. The latter therefore scales as Q −4 instead of Q −3 and the cross-section (at t = 0) correspondingly scales as We can check this analytical estimate numerically and the results are shown in Fig. 4. We plot there the coherent ρ and J/Ψ production cross-section at t = 0 scaled by Q 8 . The proposed scaling from our simple argument is not as accurate as in the case of longitudinal polarization. This is because the contribution for transversely polarized photons from large dipoles is not suppressed by large Q 2 , a consequence of the strong dependence of the overlap wavefunction in this case on the z → 0, 1 limits. Even if the Q 8 scaling is not apparent in the Q 2 range accessible at the EIC, the large suppression of the exclusive vector-meson cross-section with Q 2 is striking and can be cleanly tested.
IV. RESULTS FOR LOW Q 2 : Q 2 < Q 2 s,A Deep in the saturation region, the saturation scale is much larger than Q 2 ; in these asymptotics, we can approximate r 2 T Q 2 s,A 1. Thus the dipole-nucleus crosssection in this limit is given by dσ A dip /d 2 b T = 2, and the diffractive scattering amplitude in Eq. (1) becomes The Q 2 and A dependence of Eq. (16) are determined by the scale of the dipole radius and how that influences the overlap of wavefunctions. We expect that the dipole cross-section becomes independent of r for r ≥ 1/Q s,A ∼ A −1/6 . The overlap of wavefunctions in Eq. (16) for longitudinally polarized virtual photons is then In order to obtain some intuition as to what happens   6: The cross-section for coherent transverse vector meson production production at t = 0. It is Q 2 independent at low Q 2 . The cross-section is scaled in A by the asymptotic analytical expectation ≈ A 4/3 . when Q 2 Q 2 s , we will assume 3 that Q m q . We emphasize that due to the small mass of the light vector mesons, the results at small Q 2 are at the edge of the applicability of our weak coupling framework.
The scattering amplitude in these asymptotics has the form We will employ the identity where c = m q /Q s,A or another small number of the order of Q 2 /Q 2 s,A , M 2 V /Q 2 s,A in the limit of low Q 2 . The constant term dominates at small c; we therefore obtain It is reasonable to anticipate that the non-perturbative scale is r 1/M V . This sets the upper limit to 1/M V , and thus c = M V /Q s .
Since the constant dominates in the above expression, the only Q dependence we are left with is the overall Q scale from the virtual photon wavefunction. Hence the exclusive vector meson cross-section goes as for Q 2 s,A > Q 2 . This is demonstrated in Fig. 5, where the coherent vector meson production cross-section at t = 0 is shown as a function of Q 2 , scaled with Q −2 . The flattening of the obtained at low Q 2 demonstrates the Q 2 scaling. However, we note that this flattening for ρ takes place only at Q 2 0.1 GeV 2 , and our perturbative computation is not robust in that region as discussed at the beginning of this section.
Further, the approximation in Eq. (16) is not justified in case of J/Ψ production in realistic kinematics as its large mass and that of the charm quark limits the contribution from large dipoles to the cross-section. Thus in realistic kinematics, the Q 2 scaling for J/Ψ production at small Q 2 is obtained by again linearizing the dipole cross-section. This effectively adds two powers of r in Eq. (18). Noting that 1 c dxx 3 K 0 (x) = const + O(c 4 ), we find exactly the same scaling at low Q 2 than in case of ρ production, with the expectation that the low-Q 2 limit is reached already at Q 2 ≈ M 2 V . In the case of transversally polarized photons, the vector meson overlap in Eq. (17) gives ∼ rεK 1 (εr), where the extra power r comes from the derivative of φ T . Approximating again ε ≈ m q , the diffractive scattering amplitude is proportional to where we wrote x = m q r, and c = m q /Q s . Thus in this case we do not expect to have any Q 2 dependence at small Q 2 : The scaling of the vector meson cross-sections at low Q 2 for transversally polarized photons is shown in Fig. 6.
We will now combine the Q 2 scaling results by parametrizing the exclusive vector meson cross-section as and extract the Q slope parameter γ. This slope as a function of Q 2 is shown in Fig. 7 for both J/Ψ and ρ. The anticipated scaling changes from Q 2 to Q −6 for longitudinally polarized photons, and from Q 0 to Q −8 in case of transversally polarized photons. This scaling is seen clearly in the case of J/Ψ mesons for Q 2 ≤ 3 GeV 2 .
For the ρ, the requirement that Q m q is not satisfied in the kinematical domain where we consider our framework to be reliable. The expected asymptotics for the small-Q scaling exponents would be reached only at Q 2 ∼ 10 −2 GeV 2 where our weak coupling results are clearly not trustworthy. Nevertheless, the large variation 4 in γ with Q 2 and the qualitatively different behavior predicted between the ρ and the J/Ψ cross-sections are smoking guns that will indicate whether the dynamics of gluon saturation is at play. These results will be further corroborated by the A dependence which we will now turn to.
The only dependence on A in Eq. (16) for the low Q 2 region comes from the d 2 b T integral, which gives and the total coherent cross-section scales like A 2/3 . As can be seen in Figs. 5 and 6, this asymptotic scaling regime is not reached in realistic kinematics. For example, if we look at longitudinal ρ production at small Q 2 in the realistic kinematics of Fig. 5, the A scaling turns out to be approximately A 1.7 instead of A 4/3 . The numerically computed A dependence in the IPsat model can be expressed as Similarly, we can parametrize the total coherent crosssection as We will extract these exponents as a function of Q 2 for large nuclei. These are shown for the differential crosssection in Fig. 8 for both J/Ψ and ρ at x P = 0.001 (the energy dependence of the scaling exponents is weak). Similarly, the exponents for the total coherent crosssection are shown in Fig. 9. The analytical A 4/3 scaling for t = 0 production at low Q 2 is not reached in the kinematical domain accessible in the future electron-ion colliders. However the change in the ρ cross-section from ∼ A 1.5 to ∼ A 2 (∼ A 0.9 to ∼ A 4/3 in case of the total cross-section) when moving from low Q 2 to the dilute region is observed. In case of J/Ψ, at low Q 2 we are quite far from the analytical estimate. This is because the large mass suppresses contributions from the saturated region, and the asymptotics of dσ dip /d 2 b T = 2 corresponding to a saturated dipole amplitude in Eq. (16) is not valid.

V. CONCLUSIONS AND OUTLOOK
We demonstrated here how saturation effects significantly modify the Q 2 and A scaling properties of the exclusive vector meson production cross-section in the cross-over from the perturbative QCD regime of large Q 2 's to the saturation regime of Q 2 ≤ Q 2 s,A . In addition to analytical estimates that are valid in asymptotic kinematics, we presented numerical results for the magnitude of these effects in the kinematics relevant for the Electron-Ion Collider.
The total cross-section for exclusive vector meson cross-section (integrated over t) at low Q 2 has an A dependence between A 2/3 and A, with the numerical result closer to the latter power law dependence. This changes to A 4/3 in the pQCD regime at Q 2 Q 2 s,A . Similarly, as one goes from Q 2 Q 2 s,A to Q 2 Q 2 s,A , from perturbative QCD to deep within the saturation regime, one observes that the exclusive longitudinal ρ meson cross-section changes its Q 2 dependence from 1/Q 6 to Q 2 (1/Q 8 to Q 0 in case of transverse polarization). The scaling relations in asymptotic kinematics are summarised in Table I. The observation of such qualitative systematics would provide strong evidence for gluon saturation.
Our calculations are performed to leading logarithmic accuracy. We note that there has been much progress Longitudinal, low Q 2 Longitudinal, high Q 2 Transverse, low Q 2 Transverse, high Q 2 dσ γ * +A→V +A /dt (t = 0) I: Analytical estimates for the scaling laws for the differential cross-section at t = 0 and for the total cross-section.
in developing the saturation picture beyond leading log accuracy [51][52][53][54][55][56][57]. In particular, exclusive vector meson production at NLO has also been computed recently [29]. Despite these advances, the NLO results are not completely robust to be applied to phenomenological studies [58,59]. It will be an important topic of future research to see if the NLO results significantly modify the results given here. Finally, we note that in addition to looking at ratios of exclusive vector meson cross-sections for different nuclei, further insight can be obtained by looking at these ratios in central e + A collisions relative to those in peripheral collisions [60].