Accessing the gluonic structure of light nuclei at the Electron Ion Collider

We show how exclusive vector meson production off light ions can be used to probe the spatial distribution of small-$x$ gluons in the deuteron and $^3$He wave functions. In particular, we demonstrate how short range repulsive nucleon-nucleon interactions affect the predicted coherent $J/\Psi$ production spectra. Fluctuations of the nucleon substructure are shown to have a significant effect on the incoherent cross section above $|t|\gtrsim 0.2\,\mathrm{GeV}^2$. By explicitly performing the JIMWLK evolution, we predict the $x$-dependence of coherent and incoherent cross sections in the EIC energy range. Besides the increase of the average size of the nucleus with decreasing $x$, both the growth of the nucleons and subnucleonic hot spots are visible in the cross sections. The decreasing length scale of color charge fluctuations with decreasing $x$ is also present, but may not be observable for $|t|<1\,\mathrm{GeV}^2$, if subnucleonic spatial fluctuations are present.


I. INTRODUCTION
Electron-proton and electron-nucleus collisions can be used to precisely probe the internal structure and dynamics of protons and nuclei. Deep inelastic scattering (DIS) measurements of electrons on protons, performed at HERA, in which the electron emits a virtual photon which scatters off the target proton, have provided a detailed picture of the internal quark and gluon structure of the proton [1,2]. These measurements have revealed that at high energies (small longitudinal momentum fraction x), the proton structure is dominated by gluons.
Recently, the authors have argued that the spatial distribution of small-x gluons in the proton fluctuates eventby-event. This is evident from studying exclusive vector meson production. In coherent scattering where the target proton remains intact, the average shape of the proton is probed. In incoherent diffraction where the target dissociates, on the other hand, one is sensitive to the amount of event-by-event fluctuations [3,4] (see also [5][6][7][8][9]). These nucleon shape fluctuations have also been suggested to have a measurable effect in heavy nuclei [10].In addition to DIS experiments, the nucleon shape can be studied in hadronic collisions. In particular, the protonlead collisions at the Large Hadron Collider (LHC) have revealed unexpected collective phenomena (for a review, see e.g. [11]). One potential source of the observed collectivity is the final state response to the initial state geometry. To verify this interpretation and to disentangle it from other sources of correlations, a good understanding of the proton geometry (and nucleon geometry in the nucleus) is required. It was shown that the flow measurements in LHC proton-lead collisions are compatible with a hydrodynamically evolving Quark Gluon Plasma (QGP), initiated with a Color Glass Condensate (CGC) initial condition, only if proton geometry fluctuations are included [12,13]. Similarly, the geometry fluctuations in the nucleons were found to be important in a global analysis of lead-lead and proton-lead flow data [14]. Also, analysis of the elastic proton-proton differential cross section, measured by the TOTEM Collaboration [15] at high energy, has revealed indications for a hot spot structure of the proton [16].
The structure of light nuclei, such as the deuteron and helium at large x, is well known from low-energy scattering experiments. However, little is known about the spatial distribution of small-x gluons in these systems. To access the fundamental information on the small x gluon structure in light nuclei, as well as to provide input to models aimed at describing deuteron-gold and helium-gold collisions at the Relativistic Heavy Ion Collider (RHIC) [17][18][19], one needs new measurements possible at a future electron ion collider (EIC) [20].
In this work we study the EIC's potential to constrain the spatial distribution of small-x gluons in light nuclei via measurements of exclusive vector meson production. In particular, we compute differential coherent and incoherent cross sections for J/Ψ production within both the IPsat model [21] and the Color Glass Condensate framework, which includes the explicit solution of the JIMWLK [22][23][24][25]

evolution equations.
This article is structure as follows. In Sec. II we discuss how the deuteron and helium structure is obtained in terms of the nucleon constituents. In Sec. III it is shown how diffractive scattering processes are calculated at high energy in the dipole picture. The required dipole-nucleus scattering amplitudes encoding the target structure are obtained as discussed in Sec. IV. The resulting energy evolution for the structure of light nuclei is illustrated in Sec. V, and predictions for the future Electron Ion Collider are shown in Sec. VI. Our conclusions and outlook are presented in Sec. VII.

II. STRUCTURE OF LIGHT IONS
The density profiles of light nuclei can be calculated rather accurately using theoretical methods to describe strongly correlated quantum systems. We will extract the nucleon positions in light (A = 2, 3) nuclei from such calculations configuration-by-configuration. In the high energy scattering processes at the EIC one is sensitive to the small-x gluon distribution, about which we have to make additional assumptions given the distribution of nucleons. Future experimental data in combination with our calculations will be able to better constrain the small-x structure of light nuclei.
In order to quantify the uncertainty in the current understanding of the deuteron structure, and to study the capabilities of the EIC for constraining the model uncertainties, we apply two different deuteron wave functions in this work. First, we use the deuteron wave function from [26,27], obtained using the Argonne v18 (AV18) two-nucleon potential, including both S and P wave contributions, referred to as Argonne in this manuscript. This ab initio calculation includes both attractive and repulsive nucleon-nucleon interactions.
For comparison, we also employ the commonly used Hulthen wave function [28] in which the distribution of the proton-neutron separation d pn is parametrized as The experimentally determined parameters used in this work are a = 0.228 fm −1 and b = 1.18 fm −1 [28]. Note that these parameters are fixed by low energy data, and there is no a priori reason why the small-x gluonic distribution should resemble this precisely. The distribution of the proton-neutron separation d pn (in 3 dimensions) obtained from these two parametrizations is shown in Fig. 1. Both the Hulthen and Argonne potential wave functions produce deuterons with roughly the same root mean square size (3.93 fm in case of the Argonne potential, 4.07 fm in case of Hulthen). The largest difference between the wave functions is that the repulsive short-range nucleon-nucleon interactions suppress small proton-neutron separations in the Argonne wave function compared to the Hulthen case. We note that recently the short range correlations in deuterons (and other light nuclei) have been studied in detail, in connection with the EMC effect, see e.g. [29][30][31].
The nucleonic structure of 3 He is obtained from the Monte Carlo calculation with AV18+UIX interaction [32]. AV18 refers to the same two-nucleon potential used to obtain the deuteron configurations discussed above. In practice, we use the same database of configurations that is used in Ref. [33], and available in the PHO-BOS Monte Carlo Glauber implementation [34]. These configurations reproduce the charge radii and form factors of Helium-3, as well as the relative separation of proton pairs.

III. DIFFRACTIVE DEEP INELASTIC SCATTERING
Diffractive vector meson production is a powerful probe of the spatial structure of nuclei, as the total momentum transfer, which is the Fourier conjugate to the impact parameter, is measurable. These processes are divided in two categories: in coherent scattering the target hadron remains in its ground state, and in incoherent diffraction the target dissociates. In either case there is no net color charge transferred to the probe, which leads to an experimentally observable rapidity gap between the produced vector meson and the target or target remnants.
In the Good-Walker picture [35], the coherent cross section is obtained by calculating the average scattering amplitude for the Fock states of the probing virtual photon that diagonalize the scattering matrix. At high energy, these states are quark-antiquark states with fixed transverse separation r scattering off a fixed target configuration. The coherent cross section can be written as [5,36] where the average refers to the average over all possible target configurations, and N can be a proton or a nucleus. The scattering amplitude A γ * N →V N (x P , Q 2 , ∆) can be written as a Fourier transform from coordinate space to momentum space [36], Here, the two-dimensional vector ∆ is the transverse momentum transfer to the target, with |∆| ≈ √ −t. The transverse momentum transfer ∆ is actually the Fourier conjugate to b−(1−z)r because of the contribution from the non-forward vector meson wave function [36,37]. The impact parameter b points to the center of the dipole from the center of the target. The virtual photon to quark-antiquark dipole splitting is described by the virtual photon wave function Ψ, and the formation of a vector meson is encoded in the vector meson wave function Ψ V . In this work we use the Boosted Gaussian wave function parametrization from Ref. [36] (see also Ref. [38]). Note, that the limited knowledge of the vector meson wave function leads to significant uncertainties in the absolute normalization of the cross section (up to ∼ 50%), see e.g. Ref. [39]. The interaction of a dipole with transverse size r and impact parameter b with the target (proton or nucleus) is described in terms of the dipole amplitude N (r, b, x P ). The target structure is probed at Bjorken-x which can be interpreted as the longitudinal momentum fraction of the proton carried by the color-neutral "pomeron" exchanged with the diffractively produced vector meson (J/Ψ in this work). Here W is the centerof-mass energy in the photon-nucleon scattering and m N the nucleon mass. If we calculate the total cross section for diffractive vector meson production and subtract the coherent contribution, we obtain the incoherent diffractive cross section. Following Ref. [8] (see also [5][6][7]9]) the incoherent cross section becomes As the incoherent cross section is a variance, it measures the amount of fluctuations in the diffractive scattering amplitude. Additionally, the coherent cross section depends on the average scattering amplitude and consequently on the average target structure. These two cross sections then make it possible to constrain the eventby-event fluctuating structure of the hadron, as shown in case of protons in Ref. [3] and with heavy nuclei in Ref. [10]. There are two phenomenological corrections to the results presented above (see Ref. [36]). First, in the dilute limit where two gluons are exchanged, one actually probes a two-gluon distribution of the target. A dominant contribution to the cross section then originates from the configuration where one of the gluons is very soft, and in this limit one can relate the result to the standard collinear factorization gluon distribution function by introducing the skewedness correction [40]. However, applicability of this correction in the saturation region is not clear. The second correction originates from the fact that usually one assumes the dipole scattering amplitude to be purely real. Both of these contributions mostly affect the overall normalization of the cross section (the t and energy dependencies are weak, see e.g. [10]). As the overall normalization has a large uncertainty originating from the poorly constrained vector meson wave function, we only add these corrections approximatively. The skewedness correction is estimated as a 40% increase to the cross section. A similar 10% real part correction is applied to the results obtained using the IPsat parametrization where the dipole amplitude is purely real.

IV. DIPOLE-TARGET SCATTERING
We consider two different descriptions for the dipoletarget interaction that allow us to obtain the dipole amplitude N (r, b, x P ). These are the IPsat parametrization, in which the geometry does not evolve in x and the Bjorken-x dependence of the saturation scale Q 2 s (or density) is parametrized by fitting to HERA data. This simple parametrization is compared with an explicit Color Glass Condensate framework calculation in which the energy evolution of the Wilson lines (and consequently the target geometry) is obtained by solving perturbative evolution equations. Summaries of these two approaches are presented below.

A. IPsat model
In the IPsat parametrization the saturation scale Q 2 s depends on the impact parameter, and the dipole-proton scattering amplitude is written as where the transverse density profile function is assumed to be Gaussian: The function F contains the DGLAP evolved gluon distribution xg: The free parameters of the model (µ 2 0 , C, B p and the initial condition for the DGLAP evolution of xg) are fixed by fitting the HERA data in Ref. [41] (see also Ref. [42]).
In Refs. [3,4], this parametrization was generalized to the case where the proton has a fluctuating substructure consisting of N q "hot spots" by replacing with Here B q is the Gaussian width of each hot spot. The locations of the hot spots are sampled from a Gaussian distribution with width B qc . These parameters are constrained as in Ref. [4] by requiring a simultaneous description of the coherent and incoherent J/Ψ photoproduction data from HERA [43] at W = 75 GeV. Unlike in Ref. [4], the sampled hot spot locations are shifted to keep the center of mass of the nucleon at the origin. This keeps the deuteron and helium sizes unchanged when the fluctuations are included, but effectively makes the nucleons smaller and consequently the used parameters deviate slightly from the ones used in Ref. [4]. In case of the IPsat parametrization, we use parameters B qc = 4.5 GeV −2 and B q = 1.0 GeV −2 .
The dipole-proton scattering amplitude N discussed above can be generalized to the dipole-nucleus case as where b i are the transverse positions of the nucleons, sampled from the nucleus wave function discussed in Sec. II. Within our framework, this procedure is equivalent to summing the density profiles of the nucleons to obtain that of the nucleus.
Several pieces of data prefer the presence of additional fluctuations of the normalization of Q 2 s [3,4,44,45]. We thus allow the (squared) saturation scale Q 2 s of the individual hot spots to fluctuate around its expectation value Q 2 s following a log-normal distribution Here we use σ = 0.65 adjusted to get a better description of the small-|t| part of the incoherent J/Ψ photoproduction off protons as measured by HERA (in Refs. [3,4] with the IPsat parametrization where we did not shift center-of-mass to origin, we used comparable value σ = 0.5 determined in Ref. [45] based on observed multiplicity fluctuations in proton-proton collisions). The sampled saturation scales are then re-scaled to keep the average Q 2 s unchanged, as the expectation value of the log-normal distribution is not Q 2 s , see the discussion in Ref. [4].

B. Color Glass Condensate
The eikonal propagation of a quark through the color field of the target at transverse coordinate x is determined by the Wilson lines V (x), which describe the color rotation of the quark state. The Wilson lines at each point in the transverse plane are obtained from the Color Glass Condensate effective theory calculation [44,[46][47][48]. The details of the computation, summarized below, are exactly the same as in Refs. [3,4].
The target's color charge densities ρ a (x) in the transverse plane are assumed to be random, and sampled from the McLerran-Venugopalan (MV) model [49], where ρ a (x − , x) = 0, and where a and b are color indices. The color charge density gµ is proportional to the saturation scale Q s (x P , x, b) as Q s = cg 2 µ, with c = 0.7 when nucleon shape fluctuations are included and and c = 0.65g 2 µ without, as in Ref. [4] (see also Ref. [50]). The saturation scale Q 2 s , which is a proxy for the nucleon density, is obtained from the IPsat parametrization of the dipole amplitude, presented previously. It is defined via the relation After the saturation scale at every point in the transverse plane, Q 2 s (b), is obtained, we can sample color charges according to (12) and solve the Yang-Mills equations to determine the Wilson lines at every transverse position [48]: Here P indicates path ordering. In order to suppress long-distance Coulomb tails, an effective mass regulator m is introduced. In general one expectsm ∼ Λ QCD , and here we usem = 0.4 GeV as constrained in Ref. [4] by comparing with the HERA J/Ψ photoproduction data. All nucleon shape and density fluctuations are included via Q 2 s (b) calculated from the IPsat model, where geometry fluctuations are included as discussed in Sec. IV A.
Here we use the parameters B qc = 4.0 GeV −2 and B q = 0.3 GeV −2 , and σ = 0.5 to determine the magnitude of the Q 2 s fluctuations, as in Ref. [3]. After the Wilson lines at the initial x 0 = 0.01 are sampled as discussed above, the evolution towards smaller x is obtained by solving the JIMWLK renormalization group equation [22][23][24][25]. Here we use exactly the same setup as in Ref. [51], which we summarize below.
The JIMWLK equation describes the rapidity evolution of a Wilson line and can be written as a stochastic Langevin equation [52] d dy where t a is an SU(3) generator in the fundamental representation. The evolution (15) can be seen as a random walk in the color space, where the random noise ξ is a local Gaussian variable with variance ξ a x,i (y)ξ b y,j (y ) = δ ab δ ij δ (2) xy δ(y − y ).
The coefficient of the noise is and the kernel reads The deterministic drift term σ a x can be eliminated from the equation following Ref. [53], which avoids the requirement to evaluate Wilson lines in the adjoint representation (as appearing in (17)) and makes the numerical solution more efficient.
The JIMWLK kernel K i x has a powerlike structure, and consequently the evolution generates long-distance Coulomb tails (similar to those regulated at the initial condition by the mass parameterm in Eq. (19)). This would result in an exponentially growing cross section with rapidity [54,55], and violate the Froissart bound [56]. To avoid this we follow Ref. [55] and introduce effective confinement scale effects by using a modified kernel Here the modified Bessel function K 1 suppresses contributions at distance scale 1/m. In this work we use m = 0.2 GeV as constrained in [51] to be compatible with the HERA structure function data when the evolution is performed at fixed coupling α s = 0.21. In Ref. [51] we showed that the diffractive cross sections are sensitive to the infrared regulator m only at small t (if the strong coupling constant is adjusted to keep the evolution speed compatible with the structure function data).
Calculations are performed on a 2-dimensional lattice with transverse spacing a = 0.01 fm and L = 13 fm in case of deuterons and L = 10 fm with helium (note that larger lattices are needed with deuterons to accurately describe the tail of the distribution shown in Fig. 1). We have checked that smaller lattice spacings do not alter the results. For more details on the JIMWLK evolution and its implementation on a lattice, the reader is referred to Ref. [51].

V. ENERGY EVOLUTION OF LIGHT NUCLEI
In this section we consider the CGC framework of Sec. IV B and illustrate the evolution of individual deuteron and 3 He configurations. Results for the deuteron including proton and neutron structure fluctuations and evolution over a few units of rapidity (note that x = x 0 e −y with x 0 = 0.01 in this work) are shown in Fig. 2. The deuteron density is characterized by the trace of the Wilson line 1 − Re Tr V /N c . The evolution first washes out the proton substructure (as already noted in Ref. [55]), and eventually the nucleons grow enough to create one large region of gluon matter. However, with a typical proton-neutron transverse separation ∼ 2 fm one needs to go to very low x in order to see the two nucleons merging. Similarly, in case of 3 He, evolution for one example configuration (with nucleon substructure) is shown in Fig. 3 where very similar effects can be seen.
To quantify how the presence of two nucleons affects the deuteron evolution in comparison to the single nucleon case, we study the deuteron density evolution as a function of rapidity y = ln x0 x . In this case we fix the deuteron transverse size to a typical value d pn = 1.5 fm and quantify the density by calculating the average dipole-deuteron scattering amplitude with a fixed size dipole with |r| = 0.2 fm. Note that this is a typical scale for the dipoles contributing to the J/Ψ production [36].
For comparison, we show the case where we consider the deuteron to consist of two independently evolved nucleons. In that case, the scattering amplitude N pn for the proton neutron system at point b (where b = 0 is the center of the deuteron) is obtained following Eq. (10): where N = N p is the dipole-nucleon scattering amplitude.
The impact parameter dependence of the dipole amplitude is shown in Fig. 4, where the solid lines refer to the deuteron and dotted lines assume independently evolving nucleons. We find that at the beginning of the evolution the two nucleons are separated enough spatially and the deuteron evolution is very close to the independent nucleon case. Only at large rapidities y 4 we start to observe small deviations. We note that the rapidity at which this deviation begins depends heavily on the deuteron size d pn , and for smaller deuterons enhanced non-linear effects in the dense region start to slow down the evolution earlier.
The rapidity dependence of the dipole amplitude at zero impact parameter is shown in Fig. 5. The evolution is identical in both cases until the proton and neutron grow so much that they start to overlap significantly at y 5, when the non-linear effects start to decrease the evolution speed slightly compared with the independent nucleon case.
To demonstrate that our analysis is insensitive to the infrared regulator m, we vary m in the JIMWLK kernel (19) and extract the same dipole amplitude evolution at zero impact parameter. The value of the strong coupling constant for each infrared regulator is constrained in [51] from fits to HERA data. Besides our standard choice m = 0.2 GeV and α s = 0.21, we use m = 0.1 GeV with α s = 0.19 and m = 0.3 GeV with α s = 0.225. As can be seen from Fig. 5, the evolution is very similar in all three cases. This was not entirely obvious, because the scenario of smaller m and smaller α s for example, leads to faster evolution of the low momentum modes, which dominate the evolution of the nucleon size. For larger dipoles the sensitivity to infrared regulators should be enhanced, but contributions from large dipoles to vector meson production is suppressed exponentially.

VI. PREDICTIONS FOR THE EIC
We consider J/Ψ production in the kinematical range accessible with a future US Electron Ion Collider, where the center-of-mass energies can reach up to √ s N N = 140 Z/A GeV [57]. This allows the reach of x P values down to x P ∼ 10 −4 . . . in electron-deuteron and electron-helium collisions.

A. Short range nucleon-nucleon correlations in the deuteron
The structure of light nuclei at low energy scales is well known, but so far it has not been possible to probe the distribution of small-x gluons in e.g. the deuteron wave function. To demonstrate that the future Electron-Ion Collider can do detailed imaging of the small-x gluon distributions, we first study vector meson production off deuterons using two different realistic wave functions to describe the proton-neutron separation in the deuteron, assuming that small-x gluons are distributed around the nucleons.
The difference between the used Hulthen and Argonne wave functions is that the short range nucleon-nucleon correlations that suppress configurations where proton and neutron are close to each other (d pn 0.5 fm) are included in the Argonne wave function, see Fig. 1. For simplicity the nucleon shape fluctuations are not included in this analysis.
The obtained J/Ψ photoproduction cross sections for coherent and incoherent production at fixed x P = 0.01 are shown in Fig. 6. Here, we use the IPsat parametrization to describe the dipole-deuteron interaction. The effect of having different average density profiles for the gluon distribution results in coherent cross sections that deviate at |t| 0.3 GeV 2 . At smaller |t|, where one is sensitive to the structure at long distance scales, the spectra are identical.
Generally, the position of the first diffractive minimum corresponds to the size of the target R as t dip ∼ 1/R 2 . Here, the root mean square separations for the proton and neutron are similar in both wave functions, with the Hulthen wave function resulting in deuterons that are slightly larger by ∼ 3%. So the difference in the observed spectra must be due to the different average shape (the t spectra is the Fourier transform of the impact parameter profile, and the density profiles differ as seen in Fig. 1), not different average size. If the deuteron size is characterized by measuring the slope B D of the coherent cross section close to t = 0 (parametrizing the spectra as ∼ e −B D |t| ), we find that the spectrum is slightly steeper when the Hulthen wave function is used (B D ≈ 28 GeV −2 with Argonne potential wave function and B D ≈ 29 GeV −2 with Hulthen), consistent with the slightly larger RMS size.
The incoherent cross sections are basically identical in the studied |t| range, with the largest difference at |t| ∼ 0.1 GeV 2 , where the result with the Argonne potential wave function is below the Hulthen result by ≈ 5%. This can be understood, as the short-range correlations effectively reduce overall density fluctuations by rejecting some of the configurations where the nucleons overlap in the transverse plane. However, the effect is numerically small and most likely not observable.
To study if the difference between the two wave functions for the deuteron is washed out by the small-x evolution, we next show predictions for J/Ψ photoproduction calculated from the CGC framework using both Hulthen and Argonne wave functions to describe the proton-neutron separation at x P = 0.01. The energy evolution for each configuration is obtained by solving the JIMWLK equation. The resulting spectra at x P = 0.01 and x P = 0.0004 are shown in Fig. 7. Here, we do not include nucleon shape fluctuations for simplicity. Similar to the case of the IPsat parametrization, the coherent cross sections obtained with different wave functions deviate above |t| 0.2 GeV 2 , and the position of the first diffractive minimum is at a smaller |t| when the Argonne wave function is used. This difference is found to remain similar after the JIMWLK evolution down to x P = 0.0004. This shows that the small x evolution retains the differences between the deuteron wave functions defined at large x within the x range accessible with a future EIC.

B. Deuteron shape and its small-x evolution
First we study the effect of nucleon shape fluctuations on diffractive J/Ψ production off deuterons. In the following we will only use the Argonne potential wave function to describe the proton-neutron separation in the deuteron. Our conclusions would be similar if the Hulthen wave function were used.
Using the IPsat model with and without nucleon shape fluctuations, we calculate the coherent and incoherent J/Ψ photoproduction cross sections with detueron tar- gets at x P = 0.01. The results are shown in Fig. 8. We find that the coherent cross sections are compatible within the numerical accuracy up to |t| 1.2 GeV 2 with and without subnucleonic fluctuations (note that introducing the additional fluctuations approximately leaves the average shape unchanged). On the other hand, the incoherent cross sections start to deviate significantly at |t| ≈ 0.25 GeV 2 , similar to the case of photon-heavy nucleus scattering analyzed in Ref. [10]. Note that at smaller |t|, where we are sensitive to fluctuations at longer length scales (scale of the proton and neutron separation, sampled from the deuteron wave function), there is basically no difference in the incoherent cross sections. Next, we study the energy dependence of the cross sections. We calculate J/Ψ production at different x P values in the EIC energy range, and the results are shown in Fig. 9 for the IPsat model, and in Fig. 10 for the CGC framework, where the JIMWLK evolution equation describes the structure evolution. The results are shown with nucleon shape fluctuations (solid lines) and without (dotted lines).
We find that the position of the first diffractive minimum moves to smaller |t| as a result of the JIMWLK evolution in the Color Glass Condensate calculation, which results in growing deuteron size as illustrated in Fig. 2. On the other hand, in case of the IPsat parametrization where there is no geometry evolution the dip location is approximately constant. Similarly, the coherent cross section at small |t| gets steeper when the JIMWLK evolution is performed and remains constant in the IPsat calculation. This is demonstrated explicitly by parametrizing the coherent cross section as dσ/dt ∼ e −B D |t| and extracting the |t| slope B D as a function of x P , shown in Fig. 11. As shown in Ref. [9], the |t| slope of the incoherent cross section B incoh is controlled by the size of the object that is fluctuating. If the cross section is parametrized as e −B incoh |t| , at moderate |t| the slope is controlled by the size of the nucleon, with B incoh ∼ 4 . . . 5 GeV −2 (note that the nucleon RMS size is set by B p = 4 GeV −2 ). At large |t|, if the substructure fluctuations are included and the nucleons consist of hot spots, the slope approaches the hot spot size and B incoh ∼ 1 . . . 2 GeV −2 (recall that our hot spot size in the IPsat parametrization is B q = 1 GeV −2 ). The |t| slopes extracted from the IPsat calculation of the incoherent cross section at two different Bjorken-x values are shown in Fig. 12. As the size of the nucleons and hot spots does not depend on x, the extracted slopes are independent of x. In the CGC calculation, on the other hand, both nucleons and hot spots grow as a result of the small-x evolution. The extracted slopes in this case are shown in Fig. 13. At moderate |t| ∼ 0.2 GeV 2 the result of the x evolution is to make the spectra steeper, because of the growth of the system, both in the case with and without nucleon shape fluctuations. For any x, at large |t| nucleon shape fluctuations start to dominate if included, similar to the IPsat calculation. Without nucleon shape fluctuations, the slopes are not constant at large |t| in contrast to the IPsat calculation. This is due to the short-range color charge fluctuations in the target that contribute to the variance of the scattering amplitude. As a result of the small-x evolution, these fluctuations start to take place at shorter distance scales, and consequently the cross section falls more slowly as a function of |t| after the evolution. Because of this, the slope as a function of |t| for smaller x crosses that for the initial x at an intermediate |t| ≈ 0.6 GeV 2 .
We note that, in case of an infinite target, in the limit |t| Q 2 s and for small dipoles, the incoherent cross section approaches a power-law, which dominates over the exponentially falling geometric component at large enough |t|. When the nucleon shape fluctuations are included, the short-range color charge fluctuations do not have a large effect as the results are very similar to the case of the IPsat calculation in the studied |t| range.

C. Evolution of the Helium structure
Similar to the case of deuterons, we study the effect of nucleon shape fluctuations in 3 He by calculating diffractive J/Ψ photoproduction cross sections off helium-3 in the EIC kinematics, first using the IPsat parametrization to describe dipole-nucleon scattering. The results are shown in Fig. 14, where we again find that the effect of nucleon shape fluctuations changes the incoherent cross section at |t| 0.2 GeV 2 , corresponding to the distance scale of the nucleon substructure used in our calculations.
To study the energy evolution of the helium structure illustrated in Fig. 3, we calculate coherent and incoherent J/Ψ production cross sections at the initial condition and at smaller x values accessible at the EIC. Similar to deuterons, the position of the first diffractive minimum moves to smaller |t| as a result of the evolution, and the coherent spectra get steeper at small |t|. Inco-herent cross sections in both cases are comparable up to |t| ∼ 0.2 GeV 2 , after which the spectra are much more steeply falling in case of no nucleon substructure fluctuations.

VII. CONCLUSIONS AND OUTLOOK
We presented predictions for exclusive J/Ψ photoproduction in the EIC kinematics with deuteron and Helium-3 targets. We showed that employing two commonly used deuteron wave functions resulted in significantly different coherent vector meson spectra. This demonstrates the sensitivity for probing the transverse distribution of small-x gluons at a future EIC. In particular, diffractive vector meson measurements at an EIC could reveal deviations in the spatial gluon distributions of light nuclei from the model assumptions in this work.
By solving the small-x JIMWLK evolution equation, we predicted the energy dependence of the coherent and incoherent cross sections in the EIC energy range. The slope of the coherent |t|-dependent cross sections was found to become steeper due to the growth of the nucleus with decreasing x. The differences between the coherent cross sections computed with different wave functions at x P = 0.01 remained similar down to x P of a few times 10 −4 .
We showed that the incoherent cross section at |t| 0.2 GeV 2 is sensitive to additional nucleon substructure fluctuations, which were previously constrained by HERA data. These results are similar to what was found for heavy nuclei [10]. The small-x evolution did not significantly modify the |t| value above which the nucleon substructure affects the incoherent cross section. However, the |t| slope of the incoherent cross section also becomes steeper with decreasing x, which indicates the growth of the fluctuating constituents, namely both nucleons (|t| ∼ 0.2 GeV 2 ) and subnucleonic hot spots (|t| 0.2 GeV 2 ).
When not considering geometric subnucleon fluctuations, the slope of the incoherent vector meson spectra at |t| 0.6 GeV 2 , becomes flatter as a result of the evolution. This is because with increasing Q s , color charge fluctuations appear on decreasing distance scales. This effect is absent in calculations without explicit color charge fluctuations, as shown in case of the IPsat model.
Besides the fundamental information on the gluonic structure of light nuclei at small x that exclusive vector meson production can provide, it also has applications for the phenomenology of high energy nucleus-nucleus collisions. To interpret deuteron-gold and helium-gold measurements at RHIC, it is important to have precise knowledge of the small-x geometry of light ions, which is an input for model calculations involving hydrodynamic simulations of QGP evolution (see [11] for a review and [13,[58][59][60][61] for more recent developments).
Systematically going to heavier nuclei will be very interesting, as clustering of nucleons and other correlations, and their effect on the small x gluonic structure could be probed by a future EIC as well (see also [62][63][64] for a discussion of systematics of saturation effects in heavier nuclei).