Gluon imaging using azimuthal correlations in diffractive scattering at the Electron-Ion Collider

We study coherent diffractive photon and vector meson production in electron-proton and electron-nucleus collisions within the Color Glass Condensate effective field theory. We show that electron-photon and electron-vector meson azimuthal angle correlations are sensitive to non-trivial spatial correlations in the gluon distribution of the target, and perform explicit calculations using spatially dependent McLerran-Venugopalan initial color charge configurations coupled to the numerical solution of small $x$ JIMWLK evolution equations. We compute the cross-section differentially in $Q^2$ and $|t|$ and find sizeable anisotropies in the electron-photon and electron-$\mathrm{J}/\psi$ azimuthal correlations ($v_{1,2} \approx 2 - 10 \%$) in electron-proton collisions for the kinematics of the future Electron-Ion Collider. In electron-gold collisions these modulations are found to be significantly smaller ($v_{1,2}<0.1 \%$). We also compute incoherent diffractive production where we find that the azimuthal correlations are sensitive to fluctuations of the gluon distribution in the target.


I. INTRODUCTION
Revealing the internal structure of protons and nuclei is one of the central motivations behind the future Electron-Ion Collider (EIC) in the US [1][2][3] and similar projects proposed at CERN [4] and in China [5].These high energy deeply inelastic scattering experiments provide access to hadron and nuclear structure at small Bjorken-x, where non-linear effects [6,7] tame the growth of gluon densities and generate an emergent semi-hard scale known as the saturation momentum Q s .The existence of this scale allows for a weak coupling description of the gluon dynamics in this high density regime of Quantum Chromodynamics (QCD) in a semi-classical effective field theory (EFT) framework known as the Color Glass Condensate (CGC) [8][9][10][11][12][13][14][15][16][17].
Proton and nuclear structure in terms of fundamental quark and gluon constituents can be described in terms of parton distribution functions (PDFs).For the proton these are well known from precise structure function measurements at the HERA electron-proton collider [18].In case of nuclei, the partonic content is not as precisely known due to the lack of nuclear DIS experiments in collider kinematics [19].
To move beyond one dimensional PDFs that describe the parton density as a function of longitudinal momentum fraction, one can define, for instance, the transverse momentum dependent (TMD) parton distributions [20][21][22][23] that also include the dependence on the intrinsic transverse momentum carried by the partons, and describe their orbital angular momentum within the hadron (see e.g.[24]).Similarly, Generalized Parton Distributions (GPDs) [25][26][27] describe the spatial distribution of partons inside the hadron or nucleus [28][29][30] (see Refs. [31][32][33] for reviews on the subject).An ultimate goal in mapping the proton and nuclear structure is to determine the Wigner distribution [34][35][36] that encodes all partonic quantum information.To access these distributions, differential measurements such as two-particle or two-jet correlations are needed.These measurements have been shown to allow access to the elliptic part of the gluon Wigner distribution at small x [37][38][39], the Weizsäcker-Williams unintegrated gluon distribution [40], multi-gluon correlations [41] and gluon saturation [42][43][44][45].
The study of azimuthal correlations between the produced photon and the scattered electron in DVCS is sensitive to the gluon structure of the target.The first theoretical studies within the dipole framework in the small-x regime have been considered in Ref. [77] featuring cos(ϕ) and cos(2ϕ) modulations; the latter being generated by the elliptic gluon Wigner distribution (see also Ref. [78] for azimuthal correlations in electron quasi-elastic scattering on a large nucleus).The authors in Ref. [77] also showed the consistency of the dipole framework with the standard collinear factorization approach.In the present work, we evaluate the azimuthal correlations in DVCS and J/ψ production in the CGC EFT.In addition, we study incoherent diffractive production and demonstrate that at moderate values of |t| the azimuthal correlations are sensitive to event-by-event fluctuations of the spatial gluon distribution in the target.
This paper is organized as follows.First, in Section II, we review the calculation of exclusive cross-sections, and express the lepton-vector particle azimuthal correlations in terms of the hadronic scattering amplitudes for the production of a vector particle with polarization λ ′ in the scattering of an unpolarized nucleus with a virtual photon with polarization λ.
In Section III, we compute these amplitudes for DVCS and exclusive J/ψ production following the CGC EFT momentum space Feynman rules.We begin by considering the off-forward γ * A → γ * A amplitude, we work in the dipole picture, in which the virtual photon fluctuates into a quark-antiquark pair, scatters eikonally with the dense gluon field of the nucleus, and then recombines into a virtual photon.We express our results in terms of convolutions of light-cone wave-functions with the impact parameter dependent dipole-target scattering amplitude.We then present the explicit expressions for DVCS in Section III B and for exclusive vector meson production in Section III C. The modifications needed to compute incoherent diffraction are discussed in Section III D.
Our numerical setup is presented in detail in Section IV.We introduce both a simple impact parameter dependent Golec-Biernat-Wüsthoff (GBW) model [79] with angular correlations, and a more realistic CGC dipole computed from the spatially dependent McLerran-Venugopalan (MV) [8][9][10] initial color charge configurations coupled to small x JIMWLK 1 evolution [11][12][13][80][81][82][83].Our numerical results with predictions for the azimuthal angle correlations at the EIC are shown in Section V.As a baseline study, we present the numerical 1 JIMWLK is an acronym for Jalilian-Marian, Iancu  1. Exclusive electroproduction of a vector particle V at small x.The produced particle can be a photon or a vector meson.
results using the GBW model in Section V A. Our numerical predictions for electron-proton and electron-gold collisions from the impact parameter dependent CGC dipole amplitude are shown in Sections V B and V C for different values of virtuality Q 2 and momentum transfer |t|.
We also present predictions for the nuclear suppression factor for DVCS and exclusive J/ψ production at t = 0.In Section V D, we study incoherent diffraction with and without substructure fluctuations.We summarize our main results in Section VI and include multiple appendices (A-E) that include supplemental details.

II. EXCLUSIVE VECTOR PARTICLE ELECTROPRODUCTION
In this section we review general aspects of the computation of the vector particle (V) electroproduction crosssection in the exclusive process e(k) + A(P A ) → e(k ′ ) + V (∆) + A(P ′ A ) , shown in Fig. 1, where A can be a hadron or nucleus target with initial and final state momenta P A and P ′ A , respectively.We begin by specifying the reference frame and defining the kinematic variables of interest (see Table I for a complete list of kinematic quantities).We will express the amplitude squared in the usual way, as a product of lepton and hadron tensors, and then decompose these in a basis spanned by the longitudinal and transverse circular polarization vectors of the virtual photon exchanged in the scattering process.The advantage of using this basis is that (besides the lepton and hadron tensors having only 6 independent components, instead of 10) there is a direct connection between the polarization changing and helicity flip amplitudes and the modulations in the azimuthal angle between the electron and vector particle.
We next evaluate the components of the lepton tensor in this basis, and use symmetry arguments to deduce the generic structure of the hadron tensor and the associated amplitude for the hadronic subprocess.Using these elements, we derive a general expression for the differential cross-section of coherent diffractive vector particle electroproduction at small x, with matrix elements to be evaluated within the CGC EFT.In particular, the crosssection is differential in the azimuthal angle between the electron and the produced vector particle.The study of modulations of the cross section in this angle is the main goal of this work.

A. Reference frame and kinematics
We work in a frame in which the virtual photon and nucleon have large − and + light-cone 2 momentum components, respectively, and have zero transverse momenta: where m N is the nucleon mass.The momentum of the produced vector particle V is: In this frame, the incoming and outgoing electron carry transverse momentum k ⊥ that satisfies the relation where y and Q 2 are the inelasticity of the collision and virtuality of the exchanged photon, respectively.
Note that requiring the form of Eqs. ( 2) and (3) does not uniquely define the reference frame, due to the residual azimuthal symmetry.To completely fix the frame, we specify the azimuthal angle for the incoming electron: To characterize the collision, we introduce two invariant quantities x P and t (see Table I).At high energies, the momentum transfer squared t, from the target to the projectile is mostly transverse: Here x P describes the longitudinal momentum fraction of the nucleon carried by the effective color neutral pomeron exchange in the infinite momentum frame (IMF): outgoing four-momentum of produced vector particle V s = (P + k) 2 nucleon-electron system center of momentum squared W 2 = (P + q) 2 nucleon-virtual photon system center of momentum squared squared virtuality of exchanged photon y = P.q P.k inelasticity: in ⃗ P = 0 frame the fraction of the electron momentum carried by the virtual photon x P = (P −P ′ ).q P.q Pomeron-x: in IMF fraction of the nucleon longitudinal momentum carried by the exchanged "pomeron" 2 We define light-cone coordinates as: where v ⊥ is the two dimensional transverse vector with components v 1 and v 2 .The magnitude of the 2D vector v ⊥ is denoted by v ⊥ .The metric contains the following non-zero entries g +− = g −+ = 1 and g ij = −δ ij for i, j = 1, 2, so that the scalar product reads a

B. Lepton and hadron tensor decomposition: virtual photon polarization basis
The amplitude for the production of a vector particle V in electron-nucleus scattering is given by where λ ′ labels the polarization of the produced vector particle.The nuclear matrix element P ′ A M V,α λ ′ P A represents the amplitude for the hadronic subprocess producing a vector particle.
In light-cone gauge3 A − = 0, the tensor structure of the photon propagator has the form: By squaring the amplitude, averaging over the spin of the incoming electron, and summing over the spin of the outgoing electron and the polarizations of the particle V , we obtain the well-known decomposition for the unpolarized amplitude squared4 : 1 2 where the lepton and hadron tensors are defined by

Polarization basis
It is convenient to express the decomposition in Eq. (11) in the basis of the polarization vectors of the virtual photon.In A − = 0 gauge (and for q ⊥ = 0) they read with two-dimensional transverse vectors Here λ = 0 refers to the longitudinal polarization, and λ = ±1 denote the two transverse circular polarization states of the virtual photon.The polarization vectors satisfy the following completeness relation: where Π µα was defined for A − = 0 gauge in Eq. (10).
We can now insert 5 Eq. ( 15) on the r.h.s of Eq. ( 11) and apply the Ward identity (q α X αβ = q β X αβ = 0) to obtain an alternative expression for the unpolarized amplitude squared as 1 2 where we define the lepton and hadron tensor in the polarization basis as6 An immediate advantage of this relation is that it reduces the number of independent components of the lepton and hadron tensors from 10 to 6, which follows from the application of the Ward identities in deriving Eq. ( 16).In the remainder of the section we will see that the lepton-hadron tensor decomposition in this basis has a clear connection to the azimuthal angle correlations between the electron and the vector particle V .

Lepton tensor in the polarization basis
The lepton tensor in the polarization basis is obtained from Eqs. ( 12) and ( 17): In our choice of reference frame and because of the choice of the circularly polarized basis, the lepton tensor acquires a phase dependence on the azimuthal angle of the electron, defined in Eq. ( 6): where the elements L λ λ are given by Note that these phase structures depend on the difference between the polarization states of the virtual photon in the amplitude and the complex conjugate amplitude.

Hadron tensor in the polarization basis
The hadron tensor in the polarization basis follows from Eqs. ( 13) and ( 18): where we define The expression in Eq. ( 23) is the amplitude for the production of a vector particle V with polarization λ ′ in the collision of a nucleus with mass number A and a virtual photon with polarization λ.
In order to obtain an explicit expression for the hadron tensor, we must compute the 9 amplitudes corresponding to the different polarization states (λ, λ ′ ) of the virtual photon and the produced vector particle.This will be the subject of the next section.
However, it is possible to deduce, without explicit computation, the phase structure of the hadron tensor X λ λ.
Requiring the azimuthal rotational invariance of the decomposition in Eq. ( 16) and noticing the phase structure of the lepton tensor in Eq. ( 20), we conclude that the hadron tensor in the polarization basis decomposition must have the form: where ϕ ∆ is the azimuthal angle of the produced vector particle.For convenience we have also factored (q − ) 2 .Furthermore, by requiring the phase structure of Eq. (24), one can deduce from Eq. ( 22) that the amplitude must have the form such that The matrix elements on the right-hand side of Eq. ( 25) are independent of ϕ ∆ .

C. Cross-section, CGC average and azimuthal correlations
We will now use the results derived in the previous subsections to obtain a general expression for the unpolarized cross-section for the electroproduction of a vector particle V , and in particular the azimuthal angle correlations between the electron and the vector particle.
Our discussion thus far has been very general in terms of the nuclear matrix elements.In the coming section, we will compute these matrix elements in the CGC EFT.In the CGC one describes the target state |P A ⟩ in terms of large x stochastic classical color sources, characterized by a charge density ρ A , which generate the small x dynamical color fields.For a given operator O, one calculates the observable quantity in the CGC EFT via a double averaging procedure: 1. Compute the quantum expectation value of the operator O for a given configuration of color sources ρ A : We will see in the next section that this amounts to using modified momentum space Feynman rules that embody the effects of strongly correlated gluons inside the target at small x.
2. Perform a classical statistical average over the different color source configurations using a gauge invariant weight functional W Y [ρ A ], at a certain rapidity scale Y .The resulting quantity is the CGC average of the observable, formally defined as [11,13] The rapidity scale Y = ln(1/x) controls the separation of the small x and large x degrees of freedom, where x is typically chosen to be the longitudinal momentum fraction probed by the process (see also the discussion in Ref. [85]); we choose x = x P .
Using the following normalization for the initial and final states of the target it is easy to identify the correspondence between the nuclear matrix elements entering the amplitude squared in Eq. ( 16) (see also Eq. ( 25)) and the CGC averaged amplitude7 : Note that this prescription [86][87][88][89] of performing the CGC averaging at the level of the amplitude is specific to the case when the target remains intact after the scattering.
With these definitions in mind, it can be shown that the differential cross-section is given by where M eA,V λ ′ Y follows from Eqs. ( 9),( 15),( 23), (25) and (30) as Above F = 2k − is the electron flux factor, dk ′ and d∆ are the phase space measures of the scattered electron, and produced vector particle, respectively, defined as Recall that in our choice of reference frame k ⊥ = k ′ ⊥ .It is useful to express the phase space measure of the electron in terms of the DIS invariants: Typically, this expression is given in terms of x Bj .Here we used d ln x P = d ln x Bj .The delta function (2π)δ(q − − ∆ − ) is a manifestation of the eikonal approximation used to describe the highenergy scattering process.
Integrating the cross-section over the overall azimuthal angle and ∆ − using the delta function, we obtain where ϕ k∆ = ϕ k − ϕ ∆ is the azimuthal angle of the produced particle V with respect to the scattered electron.
Using the results for the lepton tensor in the polarization basis derived in Section II B 2, we finally obtain the differential cross-section for exclusive vector particle production correlated with the scattered electron (or electron plane): Here α em is the QED fine structure constant and 'Re' stands for the real part of the products of the CGC averaged amplitudes in the second and third lines.The explicit expressions for ⟨ M γ * A,V λ,λ ′ ⟩ Y will be obtained using the CGC EFT in Section III.
The expression for the differential cross-section above has been previously obtained at small x for DVCS in [37], and a similar expression for vector meson production in terms of its spin density matrix and helicity amplitudes has been known [90].

III. DVCS AND J/ψ PRODUCTION AT SMALL
x IN THE CGC In this section we compute the amplitudes M γ * A,V λ,λ ′ for DVCS and for J/ψ production in virtual photon-nucleus collisions at small x within the CGC EFT.
It is convenient to visualize the electroproduction of a particle V at small x in the dipole picture, in which the virtual photon γ * fluctuates into a color neutral quark−antiquark dipole q q, which subsequently scatters eikonally with the dense gluon field of the nucleus, and finally recombines into the observed final state particle V = γ, J/ψ (see Fig. 1).In the CGC EFT, the eikonal interactions (parametrically of O(1) in the QCD coupling) with the strong background classical color field are resummed into a "dressed" quark propagator, which is diagrammatically represented in Fig. 2, and can be written in momentum space as where the free Feynman propagator for a quark with fla-vor f and momentum p is given by Here m f is the mass of the quark, and i, j denote the color indices in the fundamental representation of SU(N c ).
The effective vertex for this interaction has the follow- ing expression in momentum space [91,92]: with light-like Wilson lines defined as where ρ A (z − , z ⊥ ) is the charge density of the large x color sources, t a 's represent generators of SU(N c ) in the fundamental representation and P − denotes the path ordering of the exponential along the − lightcone direction.The superscript sign(p − ) in Eq. (39) denotes matrix exponentiation: , where the latter follows from the unitarity of V (z ⊥ ).
The eikonal nature of the interaction is manifest in the conservation of the longitudinal momentum p − and the diagonal nature of the effective vertex in transverse coordinate space.
We start by computing the off-forward virtual photonnucleus amplitude M γ * A,γ * λ,λ ′ (see Fig. 3) and note that: 1. We recover the DVCS process by taking the virtuality of the final state photon equal to zero.
2. We obtain the results for J/ψ production by substituting the light-cone wave-function of the final state photon by that of the J/ψ (following the procedure in [93]).
A. The off-forward γ * A → γ * A Amplitude Let the 4-momentum vector of the outgoing virtual photon be Leading order (in the QCD coupling αs) Feynman diagram for the off-forward scattering of a virtual photon with a nucleus.The crossed circles indicate the dressed quark and antiquark propagators as defined in Eq. ( 37) with the effective vertex given by Eq. ( 39).
where we define the Recall that in our choice of frame, the virtual photon has zero transverse momentum q ⊥ = 0, while the outgoing virtual photon has nonzero transverse momentum ∆ ⊥ .
The polarization vectors for the outgoing virtual photon are then: The momentum space expression for the virtual photon-nucleus scattering matrix at small x can be written as where λ and λ ′ are the polarizations of the initial and final state photon.The trace Tr is performed both over spinor and color indices, and q f denotes the fractional charge of the quark in the loop.In writing Eq. ( 43) we introduced the short-hands: The amplitude is obtained by subtracting the non-scattering contribution (ρ a A = 0) and factoring out an overall longitudinal momentum conserving delta function: The amplitude can be written as follows: where we introduce the dipole and impact parameter vectors: In writing Eq. ( 46) we used the short-hands: The impact parameter dependent Wilson line operator in Eq. ( 46) is given by and the color-independent sub-amplitude reads with Dirac structure Performing the l − integration with the help of the delta function, and the l + and l ′+ via contour integration (see Appendix D), we obtain the simpler expression: where we introduce the light-cone longitudinal momentum fraction z = l ′− /q − , and the denominators: where we define z = 1 − z.Furthermore, following [93] we adopted the normalization8 : The result of the transverse integration over l ⊥ and l ′ ⊥ will depend on the Dirac structure A λ,λ ′ .
Let us consider the case λ = +1, λ ′ = −1, where the evaluation of the Dirac structure results in (see Appendix C): Inserting Eq. ( 55) into Eq.( 52), and evaluating the corresponding transverse integrals (See Appendix D) we find the sub-amplitude where Finally inserting Eq. ( 57) into Eq.( 46) and factoring out the corresponding phase in ϕ ∆ we obtain the amplitude where ϕ r∆ = ϕ r − ϕ ∆ .Note that this amplitude has the phase factor structure deduced in Section II B 3.
The calculations for the remaining amplitudes corre-sponding to the other polarizations states follow a similar pattern.Factoring out the corresponding phases in ϕ ∆ and a factor of q − , we find explicit expressions for , introduced in Eq. ( 25): Helicity preserving amplitudes9 : Polarization changing amplitudes (T/L → L/T): Helicity flip amplitude (T± → T∓): where we define the short-hands ζ = z 2 +z 2 , ξ = z−z, and introduce the off-forward transverse vector To obtain the matrix elements ⟨ M γ * A,γ * λ,λ ′ ⟩ Y needed to compute the differential cross-section in Eq. ( 36), we need to apply the CGG average (introduced in Eq. ( 28)) to the expressions above.This amounts to replacing the dipole operator D Y (r ⊥ , b ⊥ ) by the dipole correlator: with the operation ⟨. ..⟩Y defined in Eq. ( 28).

Accessing correlations: kinematic vs intrinsic
Eqs. ( 59)-( 63) contain phases e i(λ−λ ′ )ϕr∆ that arise from the mismatch in the polarization λ of the incoming and the polarization λ ′ of the outgoing photon states.These phases pick up different modes in the angular correlations between the dipole vector r ⊥ and the transverse momentum ∆ ⊥ .There are two sources for these correlations: 1. the intrinsic correlations between the dipole vector r ⊥ and the impact parameter vector b ⊥ in the correlator D Y (r ⊥ , b ⊥ ), combined with the Fourier phase e −i∆ ⊥ •b ⊥ , 2. the off-forward phase e −iδ ⊥ •r ⊥ .
We refer to the contribution originating from the angular correlations in D Y (r ⊥ , b ⊥ ) as intrinsic, because they are generated as a result of the target's small x structure.
Due to the even symmetry of D Y (r ⊥ , b ⊥ ) under the exchange r ⊥ → −r ⊥ , the dipole has only even harmonics10 : The intrinsic angular correlations start only at the second angular mode, and their contribution is dominant for the helicity flip term in Eq. ( 63).
In the absence of intrinsic r ⊥ , b ⊥ correlations, it is possible to have non-zero polarization changing and helicity flip amplitudes due to the presence of the off-forward phase.In this case, it is possible to show that in the dilute limit the polarization changing amplitudes are suppressed by (∆ ⊥ /Q), and the helicity flip amplitude by (∆ ⊥ /Q) 2 relative to the helicity preserving amplitudes.These scalings can also be obtained by considering the overlap between the photon and produced vector particle polarization vectors, noticing that the propagation axis for the produced particle is different from the direction of the incoming photon; therefore, we call these correlations kinematical.
We note that in the absence of the off forward transverse vector δ ⊥ , the polarization changing amplitudes vanish due to the anti-symmetry of the integrand around z = 0.5, unlike the helicity preserving and helicity flip amplitudes which are symmetric around z =0.5.

Connection to gluon GPDs
In Ref. [77] it was shown that at small x, the isotropic F 0 Y and elliptic F ϵ Y modes of the dipole correlator in momentum space are related to the unpolarized gluon GPD xH g , and the gluon transversity GPD xE T g through where Y = ln(1/x), and M is a normalization mass scale.More generally, the Fourier transform of the dipole correlator F Y can also be related to the gluon Wigner distribution at small x [37].
In the collinear limit the helicity preserving and the helicity flip amplitudes of DVCS are proportional to xH g and xE T g , respectively, as shown in Ref. [77]; thus, offering the possibility to probe gluon GPDs in exclusive particle production at small x.We emphasize that the intrinsic contribution discussed above results in a non-zero elliptic mode F ϵ Y , and consequently a non-zero transversity xE T g .
In this work we do not take the collinear limit, and instead calculate the cross-section in general (small x) kinematics.Consequently, more general structures than the two GPDs mentioned above enter in the scattering amplitudes.

B. Deeply Virtual Compton Scattering Amplitude
The amplitudes for DVCS are easily obtained by setting the virtuality of the final state photon Q ′2 = 0 in Eqs. ( 59) through (63).

Note that this immediately sets the amplitudes
±1,0 = 0, since the real final state photon can not be longitudinally polarized.The non-trivial amplitudes are given below.

Helicity preserving amplitudes:
Polarization changing amplitudes (T/L → L/T): Helicity-flip amplitudes (T± → T∓): The results above are consistent with those obtained in Ref. [77] in the massless quark limit (m f → 0), and when the dipole contains only isotropic and elliptic modulations.
We will include contributions from the light quarks u, d, s and use m f = 0.14 GeV when calculating the DVCS cross-section.
The exclusive photon production cross-section consists of the DVCS process computed above, the process of elastic electron-nucleus scattering followed by photon emission off the electron (Bethe-Heitler (BH) process), and the interference terms.In this work we only consider the DVCS process.Experimentally it is possible to at least partially distinguish the different contributions.For instance, the interference term is charge odd and can be suppressed by performing the same measurement with both electron and positron beams if sufficient precision is achieved in the experiment.The DVCS process dominates at small x P , but at the EIC energy range the BH contribution may still be non-negligible.At much smaller x P value accessible at the LHeC/FCC-he [4], the DVCS contribution should be more dominant.We also note that the two processes have a different dependence on the inelasticity y, which can further be used to enhance the DVCS signal.For a more detailed discussion related to separation of DVCS and BH contributions at the EIC, the reader is referred to Ref. [97].

C. Exclusive Vector Meson production
Following the discussion11 in Ref. [93], the amplitudes for vector meson production can be obtained by the fol-lowing replacements in Eqs. ( 59)-( 61): for the outgoing longitudinal polarization, where M V is the mass of the vector meson, and δ is a dimensionless constant (not to be confused with the off-forward vector δ ⊥ ).
We also need to perform the following replacements in Eqs. ( 60), ( 62) and ( 63) for outgoing transverse polarization.Here, we have assumed that the vector meson has the same helicity structure as the virtual photon, and introduced scalar functions ϕ T,L (r ⊥ , z) that describe the vector meson structure.
With these changes we obtain the following amplitudes: Helicity preserving amplitudes: Polarization changing amplitudes (T/L → L/T): Helicity-flip amplitudes (T± → T∓): The precise form of ϕ L,T depend on the vector meson in the final state and on the model used for the vector meson structure.In this work, we follow Ref. [93] and use the Boosted Gaussian ansatz The parameters R and N L,T in the above equation and δ in Eq. ( 73), in case of J/ψ production are determined in Ref. [93], and we use m f = m c = 1.4 GeV.We note that the use of phenomenological vector meson wave functions introduces some model uncertainty in the calculation.In addition to the phenomenological parametrizations such as the Boosted Gaussian model used in this work, there are also other approaches to describe the vector meson wave function.For example, in Ref. [98] the heavy meson wave function is constructed as an expansion in quark velocity based on Non Relativistic QCD (NRQCD) matrix elements without the need to assume an identical helicity structure with the photon.Another recent approach is to solve the meson structure by first constructing a light front Hamiltonian consisting of a one gluon exchange interaction and an effective confining potential [99,100].As shown in Ref. [98], results for the J/ψ production calculated using the phenomenological Boosted Gaussian parametrization are close to those obtained using these more systematic approaches.We note that this parametrization includes both S and D wave modes in the rest frame wave function; however, the D wave contribution results in a negligible contribution.
Exclusively produced J/ψ is experimentally measured through the dilepton channel.Similarly to the Bethe-Heitler contribution discussed in case of DVCS above, there are also QED processes with exactly the same final state (for example, elastic ep scattering, and a subsequent emission of virtual photon decaying into two leptons).We do not include this contribution in this work, but it may affect a precise extraction of the azimuthal modulations in exclusive J/ψ production.

D. Incoherent diffraction
The cross-section obtained in Eq. ( 36) corresponds to the processes where the target proton or nucleus remains in the same quantum state, and the cross-section is sensitive to the amplitudes , averaged over the target configurations.The second class of diffractive events, where the target nucleus goes out in a different quantum state A * than the initial state, is referred to as incoherent diffraction.The cross-section for eA → eA * V scattering can be obtained by first calculating the total diffractive cross-section, i.e., averaging over the target configurations at the cross-section level, and then subtracting the coherent contribution [101][102][103][104]. Consequently, the cross-section will depend on the covariances of amplitudes and thus it is proportional to where the double dipole correlator is defined as Unlike the dipole correlator in Eq. ( 64), there is no straightforward connection of the double dipole to wellknown objects in the collinear framework of QCD.Yet, it provides interesting information, as it is sensitive to the small x gluon fluctuations and event-by-event fluctuations in the color density geometry inside the nuclear target [41,65,66,104,105].See also Ref. [71] for a review.

IV. NUMERICAL SETUP
The necessary ingredient needed to evaluate the DVCS and exclusive J/ψ production cross-sections is the dipole scattering amplitude D Y (r ⊥ , b ⊥ ) introduced in Eq. (64).We perform two separate calculations to evaluate the dipole amplitude, the first using a simple GBW parametrization, the second a more realistic CGC setup.The CGC will also allow us to compute the double dipole operator that appears in the incoherent cross-sections.

A. GBW
To study how the angular modulations in the dipole amplitude affect the observable angular correlations, we use a GBW model [79] inspired parametrization as in Refs.[39,106], where angular modulation can be controlled by hand: Here ϕ rb is the angle between the impact parameter vector b ⊥ and the dipole vector r ⊥ .We set c = 0.5 to include a certain amount of angular dependence, or c = 0 when we compare with the standard GBW dipole parametrization with no dependence on the dipole orientation.The saturation scale , and transverse size of the target proton are controlled by (implicitly Y dependent) Q 2 s0 and B, respectively.As we use the GBW parametrization only to study the effect of azimuthal modulations, we do not attempt to fit the parameters precisely to HERA data, but use The cos(2ϕ rb ) modulation in (85) (in case of c ̸ = 0) results in a non-zero gluon transversity GPD, Eq.( 68), and consequently also an intrinsic contribution to the helicity flip amplitude M γ * A,V ±1,∓1 as discussed in Section III A. A particular disadvantage of this ad hoc parametrization is that very large dipoles (compared to the target size) can have unrealistic modulations e.g. when both of the quarks miss the target.Additionally, the modulation does not vanish in the homogeneous regime of small r ⊥ or b ⊥ .For a more sophisticated analytical treatment of the angular dependence in the dipole scattering amplitude, see Refs.[44,107].

B. CGC
Our main results are obtained using the CGC EFT framework that describes QCD at high energies.A particular advantage of the CGC EFT is that it is possible to calculate how the dipole scattering amplitude depends on the dipole orientation (angle ϕ rb ), and as such the effect of the elliptic gluon distribution is a prediction.For a calculation of the elliptic gluon Wigner distribution from the CGC framework, see Ref. [39].Additionally, it is possible to calculate the energy (x P ) dependence perturbatively.
The CGC framework used in this work is similar to the framework used in Ref. [63] (see also Ref. [108]), and is summarized here.The target structure is described in terms of fundamental Wilson lines V (x ⊥ ) (see Eq. ( 40)), that describe the eikonal propagation of a quark in the strong color fields of the target.As discussed in Section III, the dipole amplitude D Y is obtained as an expectation value of the correlator of Wilson lines, see Eqs. (49) and (64).Note that when evaluating the incoherent cross-sections, expectation values for the double dipole operator D (2,2) Y from Eq. ( 84) are also needed as discussed in Section III D (see also the discussion in Ref. [41]).
The Wilson lines at the initial x P = 0.01 are obtained from the MV model [8], in which one assumes that the color charge density of large x sources ρ A , is a local random Gaussian variable with expectation value zero and has the following expression for the 2-point correlator (see also Ref. [109] for recent developments) Here a, b are color indices in the adjoint representation.Following Ref. [108], the local color charge den-sity µ 2 = dx − λ A (x − ) is assumed to be proportional to the local saturation scale Q s , extracted from the IPsat parametrization [110] fitted to HERA structure function data [111,112].We note that the IPSat model, and with that our initial condition at x ≈ 0.01, includes leading order DGLAP evolution for gluons.We use Q s = 0.8g 2 µ, where the coefficient is matched in order to reproduce the normalization of the HERA J/ψ production data at x P ≈ 10 −3 .(see also [113] for a more detailed discussion relating the color charge density and the saturation scale) Given the color charge density ρ a A , one obtains the Wilson lines V (x ⊥ ) as a solution to the Yang-Mills equations as shown in Eq. (40).In the numerical implementation we include an infrared regulator m to screen gluons at large distance, and write the Wilson lines as We use m = 0.4 GeV for the infrared regulator as determined in Ref. [65].
In the IPsat parametrization, the local saturation scale ⟨Q 2 s ⟩ is proportional to the nucleon transverse density profile T p .When nucleon shape fluctuations are not included, a Gaussian profile is used where B p parametrizes the proton size.When shape fluctuations are included following Refs.[65,66], the density is written as where the proton density profile is a sum of N q = 3 hot spots whose size is controlled by B q .The hot spot positions b i ⊥ are sampled from a Gaussian distribution with width B qc .After sampling, the hot spot positions are shifted such that the center-of-mass is at the origin.Additionally, the density of each hot spot is allowed to fluctuate following the prescription presented in Ref. [65] (based on Ref. [114]) as follows.The local saturation scale of each hot spot Q 2 s is obtained by scaling ⟨Q 2 s ⟩ for each hot spot independently by a factor sampled from the log-normal distribution The sampled values are scaled in order to keep the average saturation scale intact (the log-normal distribution has an expectation value e σ 2 /2 ).Heavy nuclei are generated by first sampling the nucleon positions from the Woods-Saxon distribution.Then, the color charge density in the nucleus is obtained by adding the color charge densities of individual nucleons at every transverse coordinate.
As originally presented in Refs.[65,66], the parameters controlling the proton and hot spot sizes can be determined by fitting the coherent and incoherent diffractive J/ψ production data from HERA (see also Ref. [115] for the description of collective dynamics observed in protonlead collisions at the LHC using the constrained fluctuating proton shape).In this work, we parametrize the proton shape at the initial x P = 0.01, requiring that after the JIMWLK evolution discussed below to smaller x P ≈ 10 −3 , a good description of the HERA data is obtained. 12When no proton shape fluctuations are included, we use B p = 3 GeV −2 , and the fluctuating proton shape is described by B q = 0.3 GeV −2 , B qc = 3.3 GeV −2 and σ = 0.6.Note that shifting the center-of-mass to the origin effectively shrinks the proton, and as such both parametrizations result, on average, approximately in protons with the same size.The resulting coherent J/ψ production spectra shown in Fig. 4 that are sensitive to the average interaction with the target are almost identical.On the other hand, the incoherent cross-section from the calculation where no substructure fluctuations are included clearly underestimates the HERA data, which suggests that significant shape fluctuations are required to produce large enough fluctuations in the scattering amplitude.
The energy evolution of the Wilson lines (and consequently that of the dipole amplitude) is obtained by solving the JIMWLK evolution equation separately for each of the sampled color charge configuration.For numerical solutions the JIMWLK equation is written in the Langevin form following Ref.[116] dV The random noise ξ is Gaussian and local in spatial coordinates, color, and rapidity with expectation value zero and covariance The coefficient of the noise in the stochastic term is where The deterministic drift term reads (96) Here U (x ⊥ ) is a Wilson line in the adjoint representation of SU(N c ) and can be obtained from Eqs. ( 40) and ( 88) by replacing t a with the adjoint generators T a .For the details of the numerical procedure, in particular, related to the elimination of the deterministic drift term and expression of the Langevin step in terms of fundamental Wilson lines only, see the discussion in Refs.[63,117].The JIMWLK kernel in Eq. ( 95) includes an infrared regulator m, first suggested in Ref. [118], which exponentially suppresses gluon emission at long distances ≳ 1/m.In this work we use m = 0.2 GeV, along with a fixed strong coupling constant α s = 0.21, as constrained in Ref. [63].
As discussed above, the free parameters of the framework are constrained by the HERA J/ψ production data at W = 75 GeV [53], which corresponds to x P ≈ 10 −3 (or approximately Y = 2.3 units of JIMWLK evolution from our initial condition; the |t| dependence of x P is neglected).The agreement of the calculated J/ψ photoproduction (Q 2 = 0) cross-sections with HERA data is demonstrated in Fig. 4, where results with and without proton shape fluctuations are shown.Note that here we consider J/ψ production in γ * + p scattering.Later we will study e + p scattering where azimuthal correlations with respect to the outgoing lepton can be accessed.
In addition to vector meson production, HERA has also measured the (coherent) DVCS cross-section [72][73][74][75].The CGC framework also provides a good description of the HERA DVCS data as demonstrated in Fig. 5.When comparing J/ψ and γ production cross-sections to the HERA data, we have included an approximately 45% skewness correction [119] from Ref. [64], which takes into account the fact that in the two-gluon exchange limit, the exchanged gluons carry very different fractions of the target longitudinal momentum, and the cross-section can be related to collinearly factorized parton distribution functions.The same constant factor is included in all cross-sections shown in this work.The energy dependence of the total cross-section and t slope has been studied in [63].

V. NUMERICAL RESULTS
In this section we present numerical results for azimuthal modulations in DVCS and J/ψ production in electron-proton and electron-nucleus scattering, calculated in the energy range accessible at the future highenergy EIC.In case of electron-proton collisions, our results are shown at √ s = 140 GeV, and in case of electrongold scattering we use √ s = 90 GeV.We present results at fixed x P , which results in the inelasticity y to depend on x P , Q 2 , and |t|.We note that the azimuthal modulations and the angle independent cross-section have a different dependence on y, as apparent from Eq. ( 36).This y-dependence alone results in a suppression of azimuthal modulations with decreasing x P , especially in the case of J/ψ production.We discuss the separation of this effect from genuine JIMWLK evolution effects in Appendix B.

A. Effect of angular modulations in the GBW dipole amplitude
As discussed in Section III A, there are two contributions to the polarization changing and helicity flip amplitudes: a kinematical contribution due to the off-forward phase e −iδ ⊥ •r ⊥ , and intrinsic contribution arising from the correlations between r ⊥ and b ⊥ in the amplitude D Y (r ⊥ , b ⊥ ).
To isolate the intrinsic contribution and demonstrate the effect of the non-zero gluon transversity GPD, we first study DVCS using the GBW dipole amplitude as presented in Section IV A (note that as we calculate the cross-section in general kinematics, more general objects than the gluon GPDs introduced in Section III A are probed).This parametrization allows us to easily include the intrinsic contribution by introducing the dependence on the dipole orientation relative to the impact parameter using the parameter c introduced in Eq. (86).
The exclusive photon production (or DVCS) crosssection computed using the GBW dipole amplitude as a function of squared momentum transfer |t| ≈ ∆ 2 ⊥ is shown in Fig. 6a with no angular dependence in the dipole amplitude (c = 0) and in Fig. 6b when a significant angular modulation is included with c = 0.5.We show separately the three different contributions: Average refers to the total cross-section averaged over ϕ k∆ , which is the first line in Eq. (36).Additionally, the coefficients of cos(ϕ k∆ ) and cos(2ϕ k∆ ) in Eq. ( 36) are shown.Note that in DVCS the final state photon is real and has only transverse polarization states, and as such M γ * A,V λ,0 = 0 for all λ.As the modulation can change sign, in the logarithmic plots we show separately the positive and negative modulations.
The results are shown at Q 2 = 5 GeV 2 , and we have checked that the qualitative features depend only weakly on virtuality.We note that as there is no heavy mass scale in DVCS, this process is only marginally perturbative at moderate Q 2 especially in the limit z → 0, 1 (the so called aligned jet configuration), where there is no exponential suppression for the photon splitting into non-perturbatively large dipoles.In the GBW model, asymptotically large dipoles scatter with probability one, even when the quarks completely miss the target.Large dipoles also have a large elliptic modulation if c ̸ = 0, and consequently our GBW model calculation likely overestimates the cross-section and the intrinsic contribution to the modulations.Later when we use a CGC proton or nucleus as a target, non-perturbatively large dipoles (larger than the size of the target) do not contribute significantly.However, it is not clear which one of these setups includes a more realistic effective description of confinement scale effects.Consequently, in DVCS at small Q 2 , there may exist a non-perturbative contribution not accurately captured in either of the dipole picture calculations presented here.For a more detailed discussion of the effective implementation of confinement scale effects in the dipole picture, see e.g.Refs.[63,112,120,121].
Without the intrinsic contribution (c = 0, shown in Fig. 6a) the kinematical contribution results in nonzero cos(ϕ k∆ ) and cos(2ϕ k∆ ) modulations, that exhibit a diffractive pattern.The cos(ϕ k∆ ) modulation is generically larger as expected, as it is mostly sensitive to the polarization changing amplitude ⟨ M γ * A,V 0,±1 ⟩ Y .On the other hand, the cos(2ϕ k∆ ) modulation probes the helicity flip amplitude ⟨ M γ * A,V ±1,∓1 ⟩ Y which in the dilute limit is suppressed by an additional factor ∆ ⊥ /Q relative to ⟨ M γ * A,V 0,±1 ⟩ Y , as discussed in Section III A. Introducing an elliptic modulation in the GBW dipole significantly increases the cos(2ϕ k∆ ) modulation of the cross-section, and renders it even larger than the cos(ϕ k∆ ) modulation at intermediate momentum transfer |t|, as shown in Fig. 6b.The cos(ϕ k∆ ) modulation is only weakly modified by the modulations in the dipole amplitude, except at large |t| ≳ 1.5 GeV 2 .This can be understood based on Eq. (36).The cos(2ϕ k∆ ) contribution is directly proportional to the transverse helicity flip amplitude ⟨ M γ * A,V ±1,∓1 ⟩ Y , which is sensitive to the cos(2ϕ rb ) modulation in the dipole amplitude.On the other hand, the cos(ϕ k∆ ) modulation is dominated by , where neither of the terms are directly sensitive to the angular modulations in the dipole.
Both modulations again become more important at higher momentum transfer where the kinematical contribution is more important and the cos(ϕ k∆ ) modulation dominates.At large |t| the angular modulations change the average cross-section and the cos ϕ k∆ modulations because of two reasons.First, the helicity flip amplitude ⟨ M γ * A,V ±1,∓1 ⟩ Y has a small contribution to both of these components.Additionally, including the angular dependence in the dipole amplitude also changes the angle averaged interaction, due to the non-linear dependence on cos(2ϕ rb ) in Eq. (85).
The cross-sections with azimuthal modulations in exclusive J/ψ electroproduction at Q 2 = 2 GeV 2 are shown in Figs.6c (no angular dependence in the dipole) and 6d (with angular dependence).The large mass of the charm quark renders the kinematical contribution small as dipoles are generically smaller, which suppresses the contribution from the off-forward phase e −iδ ⊥ •r ⊥ .Consequently the relative importance of the cos ϕ k∆ modulation is significantly smaller than in case of DVCS studied above.The cos ϕ k∆ modulation is now positive at small |t|, in contrast to the DVCS case.The reason for this difference in sign is that in J/ψ production the final state can be longitudinally polarized, and terms ±1,0 ⟩ Y , that are absent in the DVCS case, give a large positive contribution to the cos ϕ k∆ modulation.
Similar to cos ϕ k∆ , the cos 2ϕ k∆ modulation is also small at small momentum transfers.If the angular modulations are included in the dipole amplitude, the cos 2ϕ k∆ modulation is significant in the |t| ≳ 1 GeV 2 region, similar to the DVCS case studied above.The smallness of the cos ϕ k∆ component compared to cos 2ϕ k∆ can potentially render the experimental extraction of the cos 2ϕ k∆ modulation more precise in J/ψ production compared to the DVCS case.
We have checked (within the GBW model) that the dependence of the azimuthal anisotropies on charm mass and the saturation scale Q s is weak for moderate values of |t|.We expect this to hold true in the full CGC computation that we outline next.

B. Proton target in the CGC
We now use a proton target described within the CGC framework, and study coherent DVCS followed by coherent J/ψ production in the EIC energy range.As discussed in Section IV B the initial condition for the JIMWLK evolution is constrained such that the coherent J/ψ production cross-section at W = 75 GeV (x P ≈ 10 −3 ) is compatible with the HERA data.Here we do not include the proton shape fluctuations as we consider the coherent cross-section, which is only sensitive to the average interaction.The incoherent cross-section, which is sensitive to fluctuations, e.g. of the proton shape, is discussed in Section V D.
As discussed in Section IV B, the dependence on the dipole orientation is calculated from the CGC framework (assuming an impact parameter dependent MV model initial condition).Consequently, the intrinsic contribution to the azimuthal correlations is a genuine prediction, unlike in the GBW model calculation presented above.Additionally, the dependence on x P is a result of the perturbative JIMWLK evolution.

DVCS
The coherent DVCS production cross-section and angular modulations as a function of squared momentum transfer |t| at x P = 0.01 are shown in Fig. 7a.The results are shown again at moderate Q 2 = 5 GeV 2 .The qualitative features of our results depend weakly on Q 2 (note that our framework is not applicable in the dilute regime of large Q 2 whose dynamics is governed by linear evolution equations, namely the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [122][123][124][125]).
The intrinsic contribution predicted from the CGC calculation results in a sizeable cos(2ϕ k∆ ) modulation in the DVCS cross-section, similar to the presented GBW model calculation with an artificial angular modulation in the dipole.The cos(ϕ k∆ ) modulation, where the intrinsic contribution is subleading, is large in the small |t| region, and suppressed strongly relative to the cos(2ϕ k∆ ) contribution at larger momentum transfers.
The extracted modulation coefficients v 1 and v 2 , defined as evolution over approximately 2.3 units of rapidity to x P = 0.001 at fixed √ s = 140 GeV in the EIC kinematics.
Both |v 1 | and v 2 are found to be sizeable, up to 5−10% in the |t| region where the coherent contribution should be more easily measurable (at very large momentum transfers the incoherent contribution will dominate and render the measurement of the coherent cross-section challenging).The elliptic modulation v 2 is suppressed approximately by a factor 2 when going from x P = 0.01 to x P = 0.001, mostly due to the JIMWLK evolution (the increase in inelasticity y with decreasing x P has a negligible effect here, as demonstrated in Appendix B).This is expected as the small x evolution results in a larger proton with smaller density gradients.In the CGC framework, these density gradients result in cos 2ϕ rb modulation in the dipole-target scattering amplitude (see Appendix E), and consequently contribute to ⟨ M γ * A,V ±1,∓1 ⟩ Y defined in Eq. ( 63), which dominates v 2 .Similarly, a small decrease in the magnitude of v 1 at small |t| is expected, as the helicity flip amplitude ⟨ M γ * A,V ±1,∓1 ⟩ Y also contributes to v 1 .As a result of the small x evolution, the dip in the cos ϕ k∆ spectrum moves to larger |t| which explains why |v 1 | increases with decreasing x P at large |t|.
Dependence on the exchanged photon virtuality is illustrated in Fig. 8, where the v 1 and v 2 coefficients are computed at the initial condition and after JIMWLK evolution as a function of Q 2 .The results are obtained by using spectra integrated over 0.2 < |t| < 0.4 GeV 2 in Eq. (97).The elliptic modulation is reduced at high virtualities.This is because the characteristic dipole size scales as r 2 ⊥ ∼ 1/Q 2 , such that large dipoles, which are most sensitive to the density gradients on the proton size scale, are suppressed at high Q 2 .These large dipoles also  97).The vn's are shown both at the initial condition (x P = 0.01, black lines) and after the JIMWLK evolution (x P = 0.001, blue lines).The bands show the statistical uncertainty of the calculation resulting from averaging over the fluctuating color charge configurations.
2.5 5.0 7.5 10.0 12.5 15.0 provide the dominant part of the intrinsic contribution to ⟨ M γ * A,V ±1,∓1 ⟩ Y .The v 1 modulation depends weakly on Q 2 , and both v 1 and v 2 are suppressed as a result of small x evolution in the studied |t| range at all Q 2 .At larger |t| the Q 2 dependence is more difficult to interpret because of the presence of a sign change in the cos ϕ k∆ modulation.At x P = 0.001, the increase in inelasticity y as a function of Q 2 dominates the virtuality dependence in the Q 2 ≳ 10 GeV 2 region.We note that the modulations vanish at the kinematical boundary where y = 1.In DVCS at x P = 0.001, this boundary is at Q 2 = 19.6GeV 2 .

Exclusive J/ψ production
We now move to the discussion of J/ψ production.Here, the heavy mass of the charm quark suppresses the kinematical contributions as already discussed in case of the GBW dipole in Section V A.
The coherent J/ψ production cross-section and different azimuthal modulation contributions are shown in Fig. 9a at x P = 0.01, corresponding to the initial condition, and in Fig. 9b at x P = 0.001 after the JIMWLK evolution, both for Q 2 = 2 GeV 2 .As a result of the growing proton size, the location of the first diffractive minimum in the angle independent coherent cross-section (and in the cos(2ϕ k∆ ) modulation) moves towards smaller |t|.As mentioned above, the heavy quark mass suppresses the kinematical contribution to the modulation, leading to a negligible cos(ϕ k∆ ) modulation.On the other hand, the cos(2ϕ k∆ ) modulation remains sizeable, of the order of a few percent, as it is dominated by the intrinsic contribution.
The v n coefficients defined in Eq. ( 97) as a function of the squared momentum transfer |t| are shown in Fig. 10.The v 2 modulation is now approximately one half of the one seen in case of DVCS.The smaller elliptic modulation compared to DVCS is expected for the same reason that larger Q 2 suppresses the modulation, as discussed above.
The v 2 modulation is again strongly suppressed with decreasing x P .Now in the J/ψ production case, the change in inelasticity y as a function of x P affects the modulation coefficients more compared to DVCS, and approximately one half of the total x P dependence seen in Fig. 10 is explained by the change in y.The other half is a result of the perturbative JIMWLK evolution resulting in a larger proton with smaller density gradients.For a detailed comparison between the JIMWLK effects and the effect of changing y, see Appendix B. The v 1 coefficient is found to be at the per mill level.The virtuality dependence of v 2 (not shown) is qualitatively similar in J/ψ production and in DVCS.Contrary to the DVCS results shown in Fig. 8, the v 1 in J/ψ production is positive and has a smaller magnitude as already observed above.At x P = 0.001, we note that the kinematically allowed region in J/ψ production is In the following, let us briefly address existing experimental data, that can be compared to our results.The H1 and ZEUS collaborations have measured [47,126] spin density matrix components in J/ψ electroproduction, that can be directly related to the v n coefficients studied here.The H1 data [47] is compatible with the s-channel helicity conservation (SCHS) assumption, in which case contributions from the polarization changing and helicity flip amplitudes are zero.For the t integrated v n coefficients in J/ψ electroproduction at low Q 2 , the H1 data corresponds to |v 2 | ≲ 5% and |v 1 | ≲ 10%, compatible with our results.
The H1 collaboration has also measured the spin density matrix elements in exclusive ρ production, with the results reported in Ref. [127].As the ρ meson is much lighter, it is sensitive to larger dipoles and especially the kinematical contribution to v n is larger.The spin density matrix elements measured by H1 correspond to tintegrated v 1 = 3 . . .15% at low and moderate Q 2 .When ρ production is calculated within our CGC framework using the wave function of Ref. [93], we find a comparable modulation.The t integrated v 2 is on the order of a few percent, which can not be accurately compared with the HERA data due to the large experimental uncertainties.Ref. [128], where this process is studied in the EIC kinematics, without considering angular correlations).The v 2 modulation is found to be extremely small due to the small density gradients, except at very large |t| > 1/R 2 A where R A is the nuclear radius.However, in this region the incoherent contribution dominates, rendering the coherent process difficult to measure.The v 1 modulation, which is dominated by the kinematical contribution, is an order of magnitude larger than v 2 at small |t|, but still at a few percent level.
The cross-section for coherent J/ψ production off gold nuclei at x P = 0.01 and Q 2 = 2 GeV 2 is shown in Fig. 11b.As a result of having a small intrinsic contribution and a large quark mass rendering also the kinematical contribution small, both the cos(ϕ k∆ ) and cos(2ϕ k∆ ) modulations are negligible, at the per mill level.In both DVCS and J/ψ production, the azimuthal modulation coefficients exhibit the same diffractive pattern as the average coherent cross-section.In contrast to e + p scattering studied previously, the signs of the modulation terms do not change at the diffractive minima.
As exclusive processes are especially powerful in resolving non-linear effects that are enhanced in heavy nuclei, we also compute the nuclear suppression factor at small |t|, defined as where V = γ, J/ψ.The nuclear suppression factors as a function of x P computed at Q 2 = 2 GeV 2 in case of J/ψ production and at Q 2 = 5 GeV 2 in DVCS are shown in Fig. 12.The results are shown in the x P range accessible at the EIC with √ s = 90 GeV.Note that the coherent cross-section scales as A 2 at t = 0 in the dilute limit [46,129], and consequently R eA = 1 in the absence of non-linear effects.As discussed in Appendix A, the polarization changing terms have a negligible contribution to the azimuthal angle independent cross-sections entering in R eA .
We predict a significant nuclear suppression in the EIC kinematics, down to R eA ≈ 0.5 in case of DVCS at the smallest reachable x P values.R eA for coherent J/ψ production is approximately 0.05 above the values for DVCS in the considered kinematics where Q 2 values are different.The obtained suppression is comparable to what was found in Ref. [112] for J/ψ production using the IPsat parametrization, where the x P dependence is parametrized and not a consequence of the perturbative small x evolution (see also Ref. [129]).The J/ψ photoproduction data from LHC suggests slightly stronger suppression R eA ∼ 0.4 . . .0.6 in this kinematical domain [130].
As small |t| also dominates the t integrated coherent cross-section, we expect non-linear effects of similar magnitude in that case too (see also discussion in Ref. [112]).The advantage of examining the ratio at t = 0 is that in that case there is no need for (x P dependent) proton and nuclear form factors used in the impulse approximation to transform the photon-proton cross-section to the photon-nucleus case in the absence of non-linear effects.This would be particularly complicated here, as the effective proton form factor includes a larger contribution from the effective projectile size in the DVCS case compared to the J/ψ production.

D. Incoherent scattering
The incoherent channel is interesting for two reasons.First, as incoherent processes are sensitive to the fluctuations of the scattering amplitude, angular modula- tions in the incoherent cross-sections can be sensitive to the details of the fluctuating target geometry.Second, the incoherent cross-section dominates at large |t| where kinematical contributions to the azimuthal modulations are more important, and as such simultaneous analysis of coherent and incoherent cross-sections can potentially be used to distinguish intrinsic and kinematical contributions.In order to access fluctuations at different length scales, we consider both J/ψ and DVCS processes.
The cross-section and its modulations for the incoherent DVCS process e + p → e + p * + γ at √ s = 140 GeV and x P = 0.01 are shown in Fig. 13.Here, the pro-ton shape fluctuations constrained by the J/ψ production data are included.Similar to coherent DVCS, especially the cos(2ϕ k∆ ) modulation is sizeable, and the cos(ϕ k∆ ) modulation changes sign at |t| ∼ 2 GeV 2 .
To study the evolution of incoherent v 1 and v 2 coefficients, and quantify the effect of proton shape fluctuations, these coefficients are shown in Figs.14a and 14b at two different x P values within the EIC kinematics, calculated with and without fluctuating proton substructure.We recall that in DVCS the change in y as a function of x P has a negligible effect on modulations as shown in Appendix B.
At small |t| ≲ 0.5 GeV 2 , where the intrinsic contribution arises from fluctuations in the angular dependence of the dipole scattering amplitude at long distance scales, the substructure fluctuations have only a small effect on the v n coefficients.At large |t| where one is sensitive to fluctuations in the amplitudes M γ * A,V λ,λ ′ at short length scales, there are significant fluctuations in the angular dependence of the dipole-target scattering amplitude when proton substructure is included, which render incoherent v 2 large as shown in Fig. 14b.Smoother density gradients obtained as a result of small x evolution suppress the dependence on the dipole orientation of fluctuations, which results in decreasing v 2 with decreasing x P .
If proton substructure is not included, the non-zero incoherent cross-section and its azimuthal modulations shown in Figs.14a and 14b are consequences of the color charge fluctuations in the target proton.These fluctuations take place at a distance scale ∼ 1/Q 2 s , which decreases with decreasing x.Consequently, since at large |t| these appearing short scale structures can be resolved by the scattering dipole, they result in larger fluctuations in the dependence on the dipole orientation and lead to increasing |v n | with decreasing x P (see also the related discussion in Ref. [62]).
Next we consider incoherent J/ψ production in e+p → e + p * + J/ψ.The cross-section and its azimuthal modulations are shown in Fig. 15a (with the spherical proton) and in Fig. 15b (with proton shape fluctuations) at the initial x P = 0.01.For comparison we again show the average coherent cross-section studied above.Although the coherent cross-section probes only the average interaction, it is also slightly altered especially at large |t| by the substructure fluctuations, as the dependence on the proton density profile given by Eq. ( 90) is non-linear.We emphasize that this dependence is an artefact resulting from our choice of the parametrization for the substructure fluctuations, and in principle it would be possible to also construct a fluctuating substructure that results in exactly the same coherent cross-section (or average dipole-target interactions).
The substructure fluctuations significantly increase both the average incoherent cross-section and the magnitude of the azimuthally dependent terms (especially the cos(2ϕ k∆ ) component) in the studied |t| range.
Additionally, the fluctuations remove the sign change from the incoherent cos(2ϕ k∆ ) modulation around |t| ≈ 1 GeV 2 , and result in an exponential incoherent spectrum instead of a power law (at |t| ≳ 1 GeV 2 ).
The extracted modulation coefficients v 1 and v 2 in incoherent J/ψ production with and without proton shape fluctuations are shown in Figs.15c and 15d.The large quark mass in J/ψ production suppresses the kinematical contribution to the polarization changing amplitudes, reducing the v 1 modulation compared to the DVCS case.Although J/ψ production is generically sensitive to shorter length scales (smaller dipoles) compared to DVCS, we again find that the modulation coefficients begin to be affected by the substructure at |t| ≈ 0.5 GeV 2 , similar to the case of DVCS.This is because the momentum transfer region where one is sensitive to the possible substructure fluctuations is mostly controlled by the size of the substructure.
For |t| ≳ 0.5 GeV 2 , we predict significantly different v 2 with and without substructure fluctuations, even with opposite signs as the substructure fluctuations remove the sign change in cos(2ϕ k∆ ) spectra.A clear effect of energy evolution is also visible, with the modulation almost disappearing in the x P range accessible at the future EIC (recall that approximately one half of the x P evolution is a result of JIMWLK dynamics; see Appendix B).This suggests that a simultaneous description of the total incoherent cross-section and its v 2 modulation could potentially be used to probe more precisely the details of the event-by-event fluctuating substructure, including its x P dependence.A more detailed analysis aimed at constraining the substructure geometry and its fluctuations is left for future work.

VI. CONCLUSIONS
We have calculated cross sections and their modulations in azimuthal angle ϕ k∆ between the exclusively produced particle and the outgoing electron in DVCS and J/ψ production in electron-proton and electron-nucleus scattering.We extracted cos(ϕ k∆ ) and cos(2ϕ k∆ ) modulations and identified two separate contributions: one of kinematical origin, and another intrinsic component, which is a result of non-trivial angular correlations in the target wave function at small x.By using dipole models with and without intrinsic contribution to the crosssection, we demonstrated that especially the cos(2ϕ k∆ ) modulation is a sensitive probe of spatial angular correlations in the gluon distribution, which in the collinear limit can be reduced to the gluon transversity GPD [77].
Using a realistic CGC based setup where the energy (or x P ) dependence of the target can be calculated, we obtained predictions for the azimuthal modulations in DVCS and J/ψ production at the future EIC.Especially for DVCS we predict a significant elliptic modulation, up to 15%, and a clear energy evolution in the EIC energy range.The modulations in J/ψ production are roughly an order of magnitude smaller, but potentially still measurable.
We have also studied incoherent DVCS and J/ψ production, and demonstrated that the momentum transfer dependence of the modulation coefficients can potentially be used to constrain the detailed properties of the target's fluctuating substructure.A more detailed analysis of this possibility is left for future work.
In the scattering of electrons with a heavy nucleus the modulations were found to be highly suppressed as a result of smaller density gradients at the probed distance scales, which strongly suppresses the intrinsic contribution to the modulation coefficients.On the other hand, the angle averaged cross sections for these processes with heavy nuclear targets are sensitive to non-linear effects, and we predict significant nuclear suppression factors as low as 0.5 in the realistic EIC kinematics.
Our results provide the first explicit predictions from the CGC EFT for the azimuthal modulations measurable at the EIC.To match the expected high precision of future EIC measurements, the presented leading order calculations (which do include a resummation of high energy logarithms α s ln 1/x, to all orders) should be promoted to next-to-leading order (NLO) accuracy.Recently, there has been tremendous progress in the field towards NLO, with the photon wave function and vector meson production impact factors calculated to NLO accuracy [98,[131][132][133][134][135].Similarly, the small x evolution equations are avail-able at this order in α s [136][137][138][139][140], and first phenomenological results have been obtained recently, e.g. for the structure functions [141].We will explore performing NLO calculations of the presented processes in future work.Additionally, it is important to study how the Bethe-Heitler contribution neglected in this work affects the azimuthal modulations in DVCS in a realistic EIC setting.It would also be interesting to compare our results to predictions based on the collinear framework of GPDs, which also incorporates DGLAP evolution. of trivial y dependence, we study in this Appendix the modulation coefficients by neglecting the JIMWLK evolution.
The elliptic modulation v 2 in DVCS and J/ψ production in electron-proton collisions at √ s = 140 GeV is shown in Fig. 17.The results at x P = 0.01 (initial condition for the JIMWLK evolution) and at x P = 0.001 (after the evolution) are identical to those shown in Section V.For comparison, we also calculate the same v 2 coefficients at x P = 0.001, but now neglecting the JIMWLK evolution.This corresponds to using the same dipole-target scattering amplitude as at the initial condition (x P = 0.01), but evaluating the process kinematics (namely the inelasticity y) at x P = 0.001.In case of DVCS, the y values are generically smaller, and as such the change in y as a function of t or x P in the studied range only slightly alters the y dependent factors in Eq. (36).If the JIMWLK evolution is neglected, v 2 at x P = 0.01 and x P = 0.001 are almost identical.Thus, nearly all the observed x P dependence is a result of the genuine JIMWLK dynamics.
In exclusive J/ψ production the large invariant mass renders y larger, and consequently the change in y with decreasing x P has a numerically larger effect.As illustrated in Fig. 17, the change in y explains approximately one half of the x P dependence seen in v 2 , which suggests that the effects of JIMWLK evolution should still be visible at the EIC.

Appendix C: Evaluation of Dirac structures
In this Appendix, we provide some useful identities for the evaluation of Dirac traces and derive Eq. (55).Some useful standard identities for gamma matrices are: Tr γ i γ j γ m γ n = 4(δ ij δ mn − ε ij ε mn ) , In deriving some of these identities we used ε jk ϵ λ,k ⊥ = iλϵ λ,j ⊥ .
We list some useful identities for the computation of the Dirac structure in Eq. (55).Recall the definition of polarization vectors in Eqs. ( 14) and (42).
For ∆ ⊥ /Q ≲ 1, J 0 (∆ ⊥ /Q) ∼ 1 and J 2k (∆ ⊥ /Q) ∼ 0 for k > 0. Thus at small to moderate values of transverse momentum of the produced vector particle, the helicity preserving and helicity flip amplitudes are most sensitive to isotropic D Y,0 (r ⊥ , b ⊥ ) and elliptic D Y,2 (r ⊥ , b ⊥ ) modes of the dipole amplitude, respectively.In this limit the polarization changing amplitude is zero.At higher momentum transfers, higher k > 2 modes of the dipole D Y,k also start to contribute.
The dipole modes in the CGC proton are shown in Fig. 18a through 18d, where we normalize the elliptic and quadrangular modes defining: with D Y,n introduced in Eq. (E1).In Fig. 18a the dipole scattering amplitude averaged over the dipole orientation ⟨D Y ⟩ is shown as a function of dipole size r ⊥ and impact parameter b ⊥ .We emphasize that when convoluted with the photon and vector meson wave functions, contributions from the largest dipole sizes are exponentially suppressed.As can be seen in Fig. 18a, this also limits the impact parameter range that gives numerically important contributions to the cross-section.
The elliptic modulation v2 at the initial condition and after the JIMWLK evolution is shown in Figs.18c  and 18d.Large enough dipoles (compared to the proton size) are required to resolve the density gradients and result in non-negligible modulations.At large impact parameters large dipoles experience strong modulations, as the scattering amplitude is large when one end of the dipole hits the center of the proton and the other one is in vacuum, and vanishes when the dipole is rotated by 90 degrees when both quarks miss the proton.However, as discussed above, contributions from this range are suppressed when convoluted with the wave functions.The modulations also vanish at b = 0 where the average dipole-target interaction does not depend on the dipole orientation.The JIMWLK evolution significantly suppresses the elliptic modulation in the studied r ⊥ , b ⊥ range, as gradients are reduced for the larger proton at smaller x [39,118].
As the odderon contribution is not included, odd harmonics vanish in our setup.The higher harmonics that contribute to the cross-section at higher momentum transfers as discussed above, are found to be small, e.g.v 4 ≲ 1% in the phenomenologically important range r ⊥ , b ⊥ < 1 fm, as shown in Fig. 18b.
FIG.1.Exclusive electroproduction of a vector particle V at small x.The produced particle can be a photon or a vector meson.
FIG.7.Different contributions to the DVCS production at Q 2 = 5 GeV 2 , and the modulation coefficients vn as defined in Eq. (97).The vn's are shown both at the initial condition (x P = 0.01, black lines) and after the JIMWLK evolution (x P = 0.001, blue lines).The bands show the statistical uncertainty of the calculation resulting from averaging over the fluctuating color charge configurations.

FIG. 9 .
FIG.9.Different contributions to the coherent J/ψ production cross-section at the initial condition x P = 0.01 and after JIMWLK evolution to x P = 0.001.Results are shown at Q 2 = 2 GeV 2 .

2 FIG. 15 .
FIG.15.Upper panel: incoherent J/ψ production cross-section and its modulations with and without proton substructure fluctuations.For comparison the average coherent cross-section is also shown.Lower panel: evolution of the vn modulation coefficients.All results are at Q 2 = 2 GeV 2 .

=+
FIG. 18.Average dipole amplitude for the CGC proton and its spatial modulations as a function of dipole size r ⊥ and impact parameter b ⊥ .

TABLE I :
Kinematic variables