Diffractive dijet production in impact parameter dependent saturation models

We study coherent diffractive dijet production in electron-hadron and electron-nucleus collisions within the dipole picture. We provide semi-analytic results for the differential cross section and elliptic anisotropy in the angle between dijet transverse momentum and hadron recoil momentum. We demonstrate the direct relation between angular moments of the dipole amplitude in coordinate space and angular moments of the diffractive dijet cross sections. To perform explicit calculations we employ two different saturation models, extended to include the target geometry. In the limit of large photon virtuality or quark masses, we find fully analytic results that allow direct insight into how the differential cross section and elliptic anisotropy depend on the saturation scale, target geometry, and kinematic variables. We further provide numerical results for more general kinematics in collisions at a future electron-ion collider, and study the effects of approaching the saturated regime on diffractive dijet observables.

We study coherent diffractive dijet production in electron-hadron and electron-nucleus collisions within the dipole picture. We provide semi-analytic results for the differential cross section and elliptic anisotropy in the angle between mean dijet transverse momentum and hadron recoil momentum. We demonstrate the direct relation between angular moments of the dipole amplitude in coordinate space and angular moments of the diffractive dijet cross sections. To perform explicit calculations we employ two different saturation models, extended to include the target geometry. In the limit of large photon virtuality or quark masses, we find fully analytic results that allow direct insight into how the differential cross section and elliptic anisotropy depend on the saturation scale, target geometry, and kinematic variables. We further provide numerical results for more general kinematics in collisions at a future electron-ion collider, and study the effects of approaching the saturated regime on diffractive dijet observables.

I. INTRODUCTION
Exclusive dijet production in coherent diffractive processes in electron-nucleus (and electron-hadron) scattering provides important information on the structure of the target in both coordinate and momentum space. In certain kinematic limits it was shown [1,2] that diffractive dijet production cross sections can be directly related to the 5-dimensional gluon Wigner distribution of the target [3,4]. This can be used to constrain both generalized parton distributions (GPDs) [5,6] and transverse momentum dependent parton distributions (TMDs) [7][8][9].
Recently, diffractive dijet production in photon-hadron collisions was studied in the dipole picture [10,11], which is an appropriate framework for high energy collisions. Interesting structures, such as diffractive dips, were found in the differential cross sections as functions of the hadron recoil momentum and the mean dijet transverse momentum. Furthermore, the dependence of the cross section on the angle between those two momenta exhibited quite a complex behavior with the sign and magnitude of the elliptic anisotropy coefficient depending on the photon polarization, virtuality of the photon, mean dijet transverse momentum, and details of the dipole model.
In this work we aim to provide analytic insight into what physical properties of the projectile and target are responsible for the observed features in the angle averaged and angular dependent cross sections, and how these features depend on the kinematic variables of the studied process. For this purpose we concentrate on relatively simple dipole models, which we extend to include a spatial dependence of the target's color charge density, and, if necessary in the model, an additional explicit correlation between the dipole orientation and the impact parameter of the collision. These features are necessary to obtain non-trivial results for the diffractive cross sections and their angular dependence.
In particular, we first study the Golec-Biernat Wusthoff (GBW) model [12,13], extended to include a spatially dependent target and angular correlations. In the limit of large photon virtuality Q 2 and/or mass of the (anti-)quark in the dijet, we evaluate the diffractive dijet production cross sections analytically, which allows for important insight into how projectile and target properties affect the experimental observables.
A somewhat more realistic description, in particular of the angular dependence, can be achieved by using an impact parameter dependent McLerran Venugopalan (IP-MV) model [14][15][16], where the angular modulation of the dipole amplitude depends on spatial gradients of the transverse color charge distribution, as previously shown in [17]. We discuss how the results are modified from the GBW model in the limit of large Q 2 and/or mass of the (anti-)quark in the dijet. However, no exact analytical result can be obtained within this model, even with the simplification of the studied limit. We thus evaluate the diffractive dijet cross sections in this model numerically, which also reveals interesting features and their dependence on the target and projectile wave function, as well as the kinematic variables. We note that more complex saturation models, such as IPSat [18,19] and b-CGC [20,21] also include impact parameter dependence, but for the sake of simplicity and being able to find analytic expressions, we stick to the two simpler models mentioned above.
Finally, by varying the saturation scale and species of the target nucleus, we analyze the effects of approaching the saturated regime on the diffractive dijet observables. We find a clear nuclear dependence on the growth of the cross section with decreasing x, as well as a characteristic dip structure of the differential cross section as a function of recoil momentum, which varies with changing saturation scale. The latter is a consequence of approaching the black disk limit, which means that the effective spatial color charge distribution of the proton deviates in-creasingly from a Gaussian with decreasing x, because of unitarity constraints.
This work should provide important guidance for future experiments at an electron ion collider (EIC) [22,23], where the spatial and momentum structure of nuclear targets can be probed to unprecedented precision.
We note that geometrical effects, similar to the ones appearing in this work, have also been studied for the case of elliptic anisotropies of gluon production in A+A [24] and p+A (p+p) [17,25] collisions. However, such processes, as well as other, like inclusive dijet production in e+p collisions [26,27] will have other sources of anisotropy, dominantly from quantum interference effects. In contrast, for the case of coherent diffractive processes as studied in this work, the geometric effects are the dominant source of anisotropy.
This paper is organized as follows: In Section II we set up the notation and review the basic ingredients of coherent diffractive dijet production in the dipole framework.
In Section III we derive a semi-analytic expression for coherent diffractive dijet differential cross sections from the dipole amplitude. In Section III A we obtain expressions for the case of a dipole amplitude without angular correlations. We then extend the analysis to dipole amplitudes with angular dependence, and derive formulas for the differential cross section and elliptic anisotropy in Section III B .
In Section IV we review the models used to specify the dipole amplitude, and in Section V we investigate analytic properties of the differential cross section and elliptic anisotropy for the two models.
In Section VI we present numerical results for the coherent diffractive dijet cross sections and their angular modulation in the IP-MV model. We present both the dependence on the dijet momentum P (Section VI A) and on the recoil momentum ∆ (Section VI B). In both cases, we study proton and gold targets. We consider different photon virtualities and saturation scales expected to be reached at a future EIC. In Section VI C we discuss observable effects of approaching the saturated regime, focusing on the differential cross section.
We summarize our results and discuss the limitations of our analysis and possible extensions of these techniques to the study of other processes in Section VII. We include multiple appendices with technical details of the calculations.

II. COHERENT DIFFRACTIVE DIJET PRODUCTION IN THE DIPOLE PICTURE
A. The dipole picture In the dipole picture, coherent diffractive dijet production in electron-nucleus (-proton) scattering can be described as a two step process: the fluctuation of the virtual photon (emitted by the electron) into a color neutral quark -anti-quark dipole and the scattering of the dipole with the target (proton or nucleus). In the Color Glass Condensate (CGC) framework at high energies (small Bjorken-x) the target is described by strong classical color fields generated by partons at larger x, which are static and localized color sources. The saturation of these fields in transverse momentum space is characterized by the saturation scale Q s which grows with decreasing x. We note that our calculations are performed at leading order. Next-to-leading order (NLO) impact factors for both diffractive dijet and vector meson production have been derived [28,29]. Together with nextto leading order corrections to the CGC evolution kernel [30][31][32][33][34][35][36], complete NLO calculations will be possible.
The fluctuation of the photon with virtuality Q 2 into a quark -anti-quark dipole of size r can be computed in quantum electrodynamics and is described by the lightcone wave functions [37][38][39][40][41]: where ε 2 f = z 0 z 1 Q 2 + m 2 f and z 0 and z 1 = 1 − z 0 are the fractions of the longitudinal photon momentum carried by the quark and anti-quark, respectively. The subscript σ i ∈ {−1, 1} refers to the helicity and β j to the color of the quark (anti-quark) (i, j ∈ {0, 1}), and m f and eZ f are the mass and electric charge of the quark of flavor f . The two dimensional vector λ denotes the transverse polarization, with λ = ±1.
The functions K 0 (z) and K 1 (z) are modified Bessel functions of the second kind, which decay rapidly as their arguments increase. For massive quarks or large photon virtualities, this implies that only small dipoles will contribute to the process to be studied, which is favorable since large dipole contributions are affected by uncontrolled infrared physics at the confinement scale [11]. Consequently, we will restrict our analysis to charm quarks.
The scattering of the dipole off the strong classical color fields is treated within the eikonal approximation, in which the transverse coordinates of the dipole are unchanged as it travels through the target, and the scattering is characterized by the color rotation of the quark and anti-quark. The color rotation of the quark via multiple scattering is encoded in the longitudinal Wilson line in the fundamental representation of SU (N c = 3) Therefore, the scattering amplitude takes a simple form in transverse coordinate space (because the dipole is initially in a color singlet state, we use the color delta function δ β0β1 to contract the inner indices of the longitudinal Wilson lines to form the product) where x 0 and x 1 are the transverse coordinates, and i and j are the color indices of the quark and anti-quark after scattering off the target, respectively. The interaction of the virtual photon (with given polarization) with the proton or nucleus is given by the scattering matrix element where I ij is the unit matrix. The non-perturbative information about the degrees of freedom of the target is encoded in the longitudinal Wilson lines.
In the next section we focus on the specific case of coherent diffractive processes, which require certain constraints on the color structure of the matrix element and the procedure for averaging over the target's color charge densities.
B. From dipole amplitude to coherent diffractive dijet cross section Before we discuss how to obtain the differential cross section from the dipole amplitude, it is important to clarify the terminology regarding our process of interest. Diffractive refers to the absence of net color exchange between the dipole and the target during the scattering; i.e., the color rotation of the quark compensates that of the anti-quark. A rapidity gap is the experimental signature for this color singlet final state, as no color string is formed between dipole and target, whose breaking would lead to particle production at intermediate rapidities. We ensure a color singlet final state by restricting to the diagonal in color space.
Coherent refers to scattering processes in which the target remains intact. Coherent processes require that the average over color charge densities is taken at the level of the amplitude. This differs from inclusive processes in which the average is taken at the level of the cross section [42][43][44].
Therefore, the object of interest in coherent diffractive dijet production is given by the color-diagonal averaged matrix element: For the sake of simplicity, we will denote the piece in front of δ ij as M(x 0 , x 1 ) T λ,L . The dipole amplitude is defined as and contains information on the spatial structure of the target.
In order to compute the cross section we express the amplitude in momentum space where p 0 and p 1 are the transverse momenta of quark jet and anti-quark jet, respectively. Before presenting the differential cross section for coherent diffractive dijet production, it is useful to introduce the set of coordinates the dipole and impact parameter vectors, respectively. Their conjugates are characterizing the mean dijet transverse momentum and momentum transfer (or nucleus recoil momentum), respectively.
In the frame in which the virtual photon has large longitudinal momentum q − , the differential cross section for coherent diffractive dijet production in electron-nucleus (-proton) scattering is then given by where dΩ = (dp − 0 /p − 0 )(dp − 1 /p − 1 )d 2 P d 2 ∆, the longitudinal momenta for quark and anti-quark jet are given by p − 0 and p − 1 , respectively. The averaged amplitude is given by Here we defined Combining above results, we arrive at expressions for the differential cross section for coherent diffractive dijet production for longitudinally and transversely polarized photons. We sum over helicities and colors (and for the transverse case over the two possible polarizations λ): where we defined ζ 2 = z 2 0 + z 2 1 . The functionsF (P , ∆) andG(P , ∆) are given bỹ These functions depend on the details of the light-cone wave functions and the dipole amplitude.

III. RELATION BETWEEN MODES OF THE DIPOLE AMPLITUDE AND THE DIFFRACTIVE DIJET CROSS SECTION
In this section we derive formulas connecting the modes of the angular correlation of the dipole amplitude, and the corresponding modes of the scattering amplitude in momentum space, and the differential dijet cross section.
If the target is isotropic, the differential cross section will only depend on the magnitudes P , ∆, and the relative angle θ P ∆ = θ P − θ ∆ . The main goal of this work is to describe the dependence of the diffractive dijet differential cross section on these three variables and provide analytic insight into the relation between the target and projectile properties and experimental observables.
In order to study the angular dependence, it is useful to decompose the differential cross section in Fourier modes The first term dσ L(T ),0 /dΩ is the differential cross section averaged over angle θ P ∆ (from now on we will refer to it simply as differential cross section). The second term dσ L(T ),2 /dΩ is the elliptic anisotropy emerging from angular correlations between P and ∆.
Any momentum correlations are encoded in the func-tionsF (P , ∆) andG(P , ∆), and they arise from the angular correlations between r and b in the dipole amplitude via Fourier transform. In order to make the connection between Fourier modes of the dipole amplitude and of the functionsF andG more explicit, we decompose the dipole amplitude into Fourier modes: We can then evaluate the modes ofF (P , ∆) and G(P , ∆) using the following mathematical relation: with equivalent expressions forG. This general result shows a one-to-one correspondence between modes of a function and modes of its Fourier transform. The derivation of Eq. (20) can be found in Appendix A.
For our particular definition ofF (Eq. (16)) one has and similar expressions for G 0 and G 2 .

A. Dipole without angular correlations
We first review results in the absence of angular correlations. In this case one hasF (P , ∆) =F 0 (P, ∆) and G(P , ∆) =G 0 (P, ∆). Therefore, using the relation in Eq. (20) for k = 0, we havẽ The formulas above combined with Eq. (14) and Eq. (15) result in where in the second expression we have used the identity for the derivative of the Bessel function of the first kind J 0 (z) = −J 1 (z). These expressions have been previously obtained in [10].

B. Dipole with angular correlations
We proceed to compute the differential cross sections in the presence of angular correlations in the dipole amplitude. This is one of the main results of our paper. We will find corrections to Eq. (23) and Eq. (24), and most importantly we will derive an expression for the elliptic anisotropy (c.f. Eq. (18)). We divide this section in two parts, discussing longitudinally and transversely polarized photons separately.

Longitudinally polarized photon
In the presence of angular correlations in the dipole amplitude, we havẽ where D 0 (r, b) and D 2 (r, b) are defined in Eq. (19). This explicitly shows that the coordinate space angular correlations in the dipole amplitude produce angular correlations in momentum space. The expressions for the differential cross section and the elliptic anisotropy for longitudinally polarized photons are then The second term in the first line constitutes a small correction to the differential cross section (averaged over angle) due to angular correlations. The second line shows the elliptic anisotropy generated by angular correlations in the dipole amplitude.
Using the expression for the gradient in polar coordinates: ∂ P =P ∂ P +θ 1 P ∂ θ , we find Keeping terms only up to the second harmonic, we have The derivatives ∂ PG0 and ∂ PG2 can be computed using the identities for the derivatives of Bessel functions: The components of the differential cross section thus take the following form As in the case of longitudinally polarized photons, the additional terms in the (angle averaged) differential cross section provide a small correction. Again, the elliptic anisotropy arises from angular correlations in the dipole amplitude.
To proceed further, one needs an explicit form of the dipole amplitude. In the next section we introduce two specific models. We will then investigate analytic properties of above cross sections in Section V. In Section VI we will evaluate the cross sections numerically and present detailed results as functions of P and ∆.

IV. REVIEW OF DIPOLE MODELS
The analytic properties of the differential cross sections and elliptic anisotropies in Eqs. (28) and (36) depend on the light-cone wave function and the dipole amplitude, which encode information on the projectile and target. Because in this work we focus on understanding the analytic structure of the diffractive dijet cross sections, we will not perform complex numerical calculations, as done e.g. in [11], but introduce relatively simple models for the dipole amplitude. We focus on dipole amplitudes of the form This parametrization is appropriate for small dipole sizes compared to the color charge density gradient of the target. The second term in the exponent contains the angular correlations between r and b, which will ultimately produce the angular correlations in the cross sections. The dipole in Eq. (37) admits a simple form for the modes D 0 and D 2 , introduced in Eq. (19). By proper projection and using the integral representations of modified Bessel functions of the first kind, I 0 (z) and I 1 (z) (Eq. (B3)), we find In the limit of small dipole sizes we have In the following we introduce explicit forms for N 0 and N 2 based on the Golec-Biernat Wusthoff and impact parameter dependent McLerran Venugopalan model.

Golec-Biernat Wusthoff model
A very simple model to describe the dipole amplitude is the Golec-Biernat Wusthoff (GBW) model [12,13], where, after introducing an impact parameter dependence, N 0 in Eq. (37) takes the phenomenologically motivated form with Q s the saturation scale at zero impact parameter, and T (b) the transverse spatial profile of the target, which we assume to be isotropic. The GBW model does not contain angular correlations between impact parameter and dipole orientation. One could add the angular correlations by hand as was done in [10] by choosing for example where c characterizes the strength of angular correlations (−1 < c < 1).

Impact parameter dependent McLerran Venugopalan model
In [17] the authors computed the dipole amplitude for an impact parameter dependent McLerran Venugopalan (IP-MV) model.
They found a dipole amplitude of the form of Eq. (37) with where m is an infrared regulator, and the factor of e in the logarithm is included to regulate the divergence for dipoles of sizes larger than 1/m. We provide details on the derivation of (42) and (43) in Appendix D.
The ellipses in Eqs. (42) and (43)  It is important to point out that the angular correlations N 2 are proportional to gradients of the color charge density. This is because the dipole will probe regions of different color charge density depending on its orientation relative to the impact parameter vector. The larger the variation of the color charge density around b, the larger the angular modulations encoded in N 2 . Furthermore, the angular correlations are suppressed for large confining scale m. As m → 0, N 2 does not diverge, but m is replaced as regulator by the finite system size, or 1/R (see details in Appendix D).

V. ANALYTIC PROPERTIES
Recently, coherent diffractive dijet production has been studied numerically. Some interesting properties of the differential cross section and elliptic anisotropy were observed [11]. Our goal is to analyze the properties of the dijet production cross sections and their angular dependence based on their analytic structure. In particular we will investigate the dependence on the saturation scale Q s , as well as the photon virtuality Q 2 , photon polarization, (anti-)quark mass m f , and the geometry of the target's color charge density.
We will begin our analysis with the GBW model for which the analytic calculations are simpler, and the results will reveal interesting properties of the differential cross section and elliptic anisotropy. We then move to the more complex IP-MV model in order to study a more realistic scenario.

A. Golec-Biernat Wusthoff model
The simplest model we consider is the GBW model. Even for this simple case the remaining integrals in Eqs. (26), (27), (34), and (35) cannot be solved analytically in general. However, in the limit where the saturation scale Q s is much smaller than the photon virtuality or quark mass, analytic expressions can be found. We present results in this limit in the next section and discuss the consequences of relaxing them thereafter.

Large photon virtuality or massive quarks: Qs ε f
If Q s ε f the Bessel functions K 0 and K 1 suppress contributions from r 1/Q s in the r-integrals in Eqs. (26), (27), (34), and (35). For r 1/ε f 1/Q s , one can expand the dipole to quadratic order (Eq. (39)) and obtain In this limit the dependence of the dipole amplitude on r and b factorizes, and one finds the following expressions for the differential cross sections (we provide further details on the calculation of the following results in and the elliptic anisotropies As expected from the factorization in b and r, in this limit the P and ∆ dependencies factorize, and both the differential cross section and elliptic anisotropy grow as Q 4 s . This growth will eventually be tamed by saturation effects as we discuss in the next subsection. The P -dependence probes the projectile -it is sensitive to the photon's polarization and virtuality, as well as the quark mass encoded in ε f . The P -dependence of the differential cross section for longitudinally polarized photons develops a dip at P = ε f . The location of the dip depends directly on the virtuality Q 2 and the mass m f of the quark. This feature is absent in the transversely polarized photon case, which has two contributions, only one of which exhibits a dip. The elliptic anisotropies for both polarizations change sign at P = ε f , a feature also observed in the IPSat model in [11]. The sign of the elliptic anisotropy in the transverse case depends on the ratio ε f ζ/m f . At large P the differential cross section and elliptic anisotropy decay with the power law 1/P 8 , while for small P the former is constant, and the latter vanishes. Similarly, at Q 2 P 2 , both cross sections scale as Q −6 , while for smaller Q 2 P 2 , they scale as Q 2 . Both elliptic anisotropies scale as Q −8 for Q 2 P 2 and Q 2 for Q 2 P 2 in the GBW model. We note that detailed scaling relations with Q 2 and the mass number A for the related process of diffractive vector meson production were discussed in [45]. Information on the target is encoded in the ∆-dependence of the differential cross section, which involves the Fourier transform of the color charge density profile. The ∆-dependence of the elliptic anisotropy provides access to target properties via the productT (∆)T 2 (∆). This particular dependence on the target geometry is a consequence of how we introduced the angular dependence in the GBW model in Eq. (41). In the more realistic IP-MV model discussed below, the elliptic anisotropy will only depend onT (∆).
2. Approaching the saturated regime: Qs ∼ ε f The case where the saturation scale is of the same order as ε f is more difficult to study analytically, as one cannot expand the dipole amplitude to quadratic order, because dipoles of size r ∼ 1/Q s will contribute to the cross sections. One should use the full expressions in Eq. (38) in which the r and b dependencies do not factorize. Therefore, the P -dependence and ∆-dependence of the differential cross section and elliptic anisotropy no longer factorize either.
For fixed dipole size r, the dipole amplitude no longer grows as Q 2 s , but it saturates to unity, thus slowing down the growth of the differential cross section and elliptic anisotropies as the saturation scale Q s is increased.
While in the high virtuality or large quark mass limit the dominant momentum scale is ε f , we now have a competition between the two scales Q s and ε f . This will be reflected in the Q s -dependence of observables. For instance, the dip in P for the differential cross section in the longitudinally polarized case and the change in sign in the elliptic anisotropy will also be sensitive to the saturation scale Q s . Their locations will shift to larger values of P as the saturation scales increases. This can also be justified mathematically by observing that the location of the maximum of the product of dipole amplitude and light-cone wave function for longitudinally polarized photons, will shift towards smaller values of r as Q s increases. Thus, the Fourier transform of this product will have a zero at a larger value of P , causing the shift of the dip. Another interesting feature is that the dipole amplitude is no longer proportional to the target color charge profile T (b), thus producing a more complex ∆dependence of the differential cross section and elliptic anisotropy. The ∆-dependence will not only depend on the geometry of the profile but also on saturation effects.

B. Impact parameter dependent McLerran
Venugopalan model In this section we consider an approximation based on a parametric estimate of the effect of the logarithm in the IP-MV model, which distinguishes the outcome of the model from the result in the previously discussed GBW model. To gain more insight into the features of this model, a numerical study is necessary, which will be presented in the next section.
1. Large photon virtuality or massive quarks: Qs ε f , As done in the GBW model above, we expand the dipole for small sizes r and arrive at The presence of the logarithmic factor in D 0 makes the analytic evaluation ofF 0 (26) and ∂ PG0 (34) difficult. In Appendix C we argue that the effect of the logarithmic factor is to enhance the value ofF 0 and ∂ PG0 and to shift the zero ofF 0 to a larger value of P (See Eqs. (C16) and (C17)). We arrive at approximate expressions for the differential cross sections: where ξ characterizes the shift, and C 1 and C 2 are the enhancement factors. In Appendix C, we show that ξ > 1 and C 2 > C 1 > 1.
Similarly, one finds that the elliptic anisotropies are given by As in the GBW model for Q s ε f , the P and ∆ dependencies factorize. The P -dependence is sensitive to the light-cone wave function and the logarithm in the dipole amplitude, as manifest in the dependence on ξ. In the longitudinal case, the differential cross section displays a dip at P = ξε f (greater value of P compared to GBW -a similar behavior was observed in the location of the dip in the IPSat vs CGC results in [11]). As in the GBW model, the transverse cross section does not display a dip.
An interesting feature of the IP-MV model in this limit is that the elliptic anisotropy is highly sensitive to the infrared regulator m. One should mention that these formulas break down when m ∆, since Eq. (43) has been obtained assuming small momentum transfer kicks to the dipole (See also last paragraph of Appendix D).
The elliptic anisotropy for the longitudinal photon changes sign at P = ξε f (greater value of P compared to GBW, the shift was also seen in [11] from IPSat to CGC.).
The behavior of the elliptic anisotropy for the transverse polarization is more subtle. Let us assume that ξ ∼ 1.7 and This implies that at low Q 2 , the elliptic anisotropy is positive for small P . In contrast, for sufficiently large Q 2 , the anisotropy at small P is negative. At large P we have For our choice of ξ and C 2 /C 1 , this condition is always satisfied, and thus one expects the elliptic anisotropy to remain positive. Therefore, for transverse polarization we expect a change in sign in the elliptic anisotropy from negative to positive for high virtuality Q 2 . For low Q 2 , we expect the elliptic anisotropy to be positive for all P . Such behavior has been observed in the recent numerical study in [11]. We observe the same scaling with Q 2 and P for the differential cross section and elliptic anisotropy as in the GBW model. The ∆-dependence of the differential cross section is similar to that of the GBW case; the Fourier transform T (∆) of the color charge density appears and introduces the sensitivity to the details of the geometry of the target's color charge density. The elliptic anisotropy differs from that of GBW, in which the angular correlations where included by hand. In the IP-MV case, where the angular correlations emerge as a result of the finite size gradients of the profile, we only find a dependence oñ T (∆), not onT 2 (∆). Eqs. (53) and (54) also show that as the momentum transfer ∆ goes to zero (exact correlation limit), the elliptic anisotropy vanishes. The position of the maximum as a function of ∆ will depend on the details of the profile.

Approaching the saturated regime Qs ∼ ε f
Similarly to the case of the GBW-model, we expect the P -dependence and ∆-dependence of the differential cross section and elliptic anisotropy not to factorize, and their growth with Q s to slow down, as we move away from the limit discussed above. Compared to that limit, the P -dependent longitudinal differential cross section will develop a dip at larger values of P , and the elliptic anisotropy will change sign at larger values of P as well.
The conditions for the elliptic anisotropy for the transversely polarized photon in Eqs. (55) and (56) will be modified as Q s increases. The numerics will show that the change in sign in this case will happen for larger virtuality Q 2 . The ∆-dependence of the differential cross section and elliptic anisotropy will encode both the geometry of the target color charge profile and the saturation scale, both interlaced by the exponent in the dipole amplitude. These effects, due to the emergence of saturation, will be studied numerically in the next section.

VI. NUMERICAL RESULTS
In this section we numerically evaluate the semianalytic formulas for the differential cross section dσ 0 /dΩ and the elliptic anisotropy dσ 2 /dΩ in Eqs. (28) and (36) using the dipole amplitude of the impact parameter dependent MV model of Eqs. (42) and (43) together with the projection formulas in Eq. (38). In the first two subsections we study the P and ∆ dependence. We perform our analysis for both protons and gold nuclei, and two different photon virtualities, Q 2 = 1 GeV 2 and 10 GeV 2 . We further employ m f = 1.28 GeV for the mass of the charm quark. We choose m = 0.4 GeV as our infrared regulator. For simplicity, our analysis is performed at fixed z 0 = z 1 = 0.5, which dominates the bulk of the cross section and could be fixed in experiments as well.
For our study of the proton we use Q 2 sp = 0.5 GeV 2 and a Gaussian target profile with R p = 0.4 fm, which is the gluonic radius of the proton. The normalization is chosen such that T p (b) = 1 at b = 0. In this case the value of Q s quoted is that in the center of the proton. For our analysis of the nucleus the thickness function T A in the transverse plane is obtained by the integration of a Woods-Saxon distribution along the longitudinal direction where the Woods-Saxon distribution is given by    and N A is chosen such that T A (0) = 1. For a gold nucleus (A = 197), we choose R A = 6.37 fm and a A = 0.535 fm. The saturation scale is Q 2 sA = 1.09 GeV 2 . Once again, the normalization is such that the Q s quoted is that in the center of the nucleus. The relation between proton saturation Q 2 sp and nuclear saturation scale Q 2 sA is discussed in Appendix E. We will discuss how to use a varying mass number and Bjorken-x to uncover effects of saturation in the differential dijet cross section in Section VI C.

A. P dependence
In this section we study the P dependence of the differential cross section and elliptic anisotropy at fixed momentum transfer ∆ = 0.1 GeV. In all plots we display the longitudinal, transverse, and total cross sections separately.

Differential cross section dσ0/dΩ
We first study the differential cross sections, shown in Figs. 1, 2, 3, and 4. The most salient feature is the dip in the longitudinal cross section in all four cases. The location of the dip shifts to larger momentum with increasing photon virtuality Q 2 . The location of the dip is at larger momentum P for the gold target, compared to the proton, because the saturation scale Q s is larger for the former, and in general both Q s and Q 2 affect the P -dependence of the cross section, as discussed in the previous section.
In all four cases the transverse components dominate for most values of P . However, we see that the difference between the longitudinal and transverse cross sections decreases with increasing photon virtuality. At small P the longitudinal component dominates for Q 2 = 10 GeV 2 The cross sections decrease with photon virtuality Q 2 as this restricts the size of the dipoles contributing to the scattering. Furthermore, they are more than two orders of magnitude larger for gold nuclei compared to protons, because of the larger size of the target and the increased saturation scale (Q 2 s ∼ A 1/3 ). Observe that despite differences in the specific details, i.e., the precise location of the dip, ordering of transverse to longitudinal components, etc., the structure of the results is similar in all four cases, owing to the fact that the P -dependence is most sensitive to the form of the lightcone wave functions (projectile) and less sensitive to the target under consideration.
The magnitude of the elliptic anisotropy follows a similar pattern to that of the differential cross section, de-creasing with increasing Q 2 , and increasing with larger target size and saturation scale Q s . In all four cases one observes a maximum at P ≈ 1.0 GeV, with only a weak dependence of this location on virtuality or saturation scale (which varies with the choice of target). As expected the elliptic anisotropy vanishes at large values of P and at P = 0 GeV.
For the longitudinally polarized photon, the elliptic anisotropy changes sign in all four cases. As anticipated from our analytic investigation, this happens at values of P that coincide with the location of the dip in dσ 0 /dΩ, which shifts to larger values of P as the virtuality Q 2 and saturation scales Q s increase.
For the transversely polarized case, at low virtuality Q 2 the elliptic anisotropy remains positive (Figs. 5 and 7), while at large virtuality Q 2 , the elliptic anisotropy changes sign (Fig. 6), as discussed in the last paragraph of Section V B. Since Q s is larger for the gold nucleus, the used Q 2 is not large enough to cause a change of sign in the result shown in Fig. 8. However, as one increases the virtuality further the anticipated sign change will ap- pear: We have checked that dσ T,2 /dΩ changes sign as a function of P for Q 2 12 GeV 2 . The transversely polarized contribution to the elliptic anisotropy dominates at low Q 2 , but at high Q 2 we observe that the longitudinal piece starts to dominate for values of P below 2 GeV.
The order of magnitude for the relative total elliptic anisotropy is about 0.1% for the chosen value of ∆ = 0.1 GeV, in agreement with the results in [11].

B. ∆ dependence
In this section we study the ∆ dependence of the differential cross section dσ 0 /dΩ and elliptic anisotropy dσ 2 /dΩ at fixed P = 1.0 GeV. In all plots we display the longitudinal, transverse, and total cross sections separately.

Differential cross section dσ0/dΩ
We present the ∆-dependence of the cross section for diffractive dijet production off a proton in Figs. 9 and 10 for Q 2 = 1 and 10 GeV 2 , respectively. In both cases one observes diffractive dips at values comparable to the inverse size of the proton 1/R p . This is not purely a geometric effect, but relies on the presence of saturation. For example, the results cannot be anticipated from the non-saturated regime (quadratic expansion of the dipole amplitude, e.g. Eqs. (51) and (52)), where one would obtain a smooth ∆-dependence from the used Gaussian profile in coordinate space. In fact, this behavior cannot be obtained in any finite order Taylor expansion of the dipole amplitude. It is a consequence of the (resummed) multiple scattering of the dipole, or unitarity of the dipole amplitude [10,46]. Consequently, the location of the dip depends on Q s , which we will demonstrate explicitly in the next subsection. As before, one further observes in Figs. 9 and 10 that with increasing Q 2 the difference between cross sections for transverse and longitudinal polarizations decreases.
Results for a gold nucleus target are presented in Figs. 11, and 12. The various diffractive dips occur at multiples of ∼ 3/R A ∼ 0.1 GeV (the factor of 3 is closely related to the zeroes of the Bessel function J 0 ), where R A is the size of the nucleus. The locations of the diffractive dips also depend on a A (the skin depth of the nucleus). Unlike for the (Gaussian) proton, these dips are present due to the geometry of the nucleus. Nonetheless, their locations are also modified by the unitarization of the dipole amplitude. We have checked that in the absence of saturation (expansion to quadratic order), the dips shift to larger values of ∆. At large Q 2 the difference between the longitudinal and transverse component of the differential cross section is reduced.

Elliptic anisotropy dσ2/dΩ
The dependence of the elliptic anisotropy on the momentum transfer ∆ is shown for a proton target in Figs. 13 and 14 for Q 2 = 1 and 10 GeV 2 , respectively. Unlike the differential cross section, here we only show results for ∆ ≤ 1.0 GeV. This is because the validity of our approximation breaks down for the elliptic anisotropy at large ∆ (see the discussion at the end of Appendix D). The anisotropy increases rapidly with ∆ and reaches a maximum at ∆ ≈ 0.5 GeV, which is approximately the inverse gluonic size of the proton. At Q 2 = 10 GeV 2 , the longitudinal component dominates the elliptic anisotropy, while at the smaller Q 2 = 1 GeV 2 the transverse component dominates. This is in agreement with Figs. 5 and 6.
We conclude our discussion with the gold nucleus, shown in Figs. 15 and 16. The behavior resembles that of the differential cross section, except that a global max-  imum appears at ∆ ∼ 0.05 GeV, which is close to the inverse size of the nucleus. An interesting feature is the presence of regions where the elliptic anisotropy is negative. This is an effect of the unitary (exponentiated scattering) of the dipole amplitude, as such regions will not be present if the the dipole was expanded to quadratic order. This can be seen directly in Eqs. (53) and (54), which show that the elliptic anisotropy is proportional to ∆ 2 |T (∆)| 2 , meaning that no sign change will occur with varying ∆.
Our results on the P -and ∆-dependence of the elliptic anisotropy presented above provide important guidance for future experiments. Note that the elliptic anisotropy varies strongly as a function of both P and ∆. Consequently, the choice of the kinematics of the dijet will affect how well the anisotropy can be measured experimentally. For example, for a proton target, choosing P ≈ 1 GeV and ∆ ≈ 0.5 GeV is predicted to maximize the elliptic anisotropy, while for a gold target a smaller ∆ ≈ 0.05 GeV would be preferred.

C. Observable effects of approaching the saturated regime
In this section we study effects of saturation that are potentially observable in measurements of diffractive dijet production at a future electron ion collider. We focus our attention on the differential cross section.
In Fig. 17 we show the x-dependence of the differential cross section at fixed values of P = 1.0 GeV and ∆ = 0.1 GeV for different targets: proton, copper (R A = 4.163 fm, a A = 0.606), and gold. We normalized the differential cross sections by their values at x 0 = 0.01. In the semi-analytic model used here, the evolution in x solely affects the value of Q s . We thus simply vary Q s and relate it to x using the parametric relation For the initial saturation scales Q 2 s0 of the proton, copper and gold nucleus we chose: 0.3, 0.43, 0.65 GeV 2 , respectively. The relation between the proton saturation scale Q 2 sp and the nuclear saturation scale Q 2 sA is given in Appendix E.
The reference dotted line shows the expected evolution in the absence of saturation. In that case, the differential cross section grows with Q 4 s (see Eqs. (51) and (52)). The results show the slowdown of the growth of the differential cross section in response to saturation effects. These set in earlier for the denser targets because of their larger saturation scale Q 2 s ∼ A 1/3 for any given x. C.f. the discussion in Sec. V B 2.
In Fig. 18 we show the differential cross section as a function of ∆ for a proton target at different values of the saturation scale Q s . The narrow lines denote a "nonsaturation" model where the dipole amplitude is not exponentiated (expanded to quadratic order). The figure shows that the differential cross sections are smaller for the case including saturation (thick lines). This effect is more pronounced for larger saturation scales as expected.
The more prominent feature is the dependence of the location of the diffractive dip on the value of Q 2 s . As the saturation scale increases, the dip shifts to lower values of ∆. This effect has been observed in [10] and explained by the fact that the effective spatial shape of the proton is non-Gaussian. This happens because the center of the proton is approaching the black disk limit (mathematically, this happens because the Gaussian thickness functions appears in the exponential of the dipole amplitude). The dip is absent in the non-saturated case, as the differential cross section here is proportional to the square of the Fourier transform of the Gaussian profile, which does not have a dip structure (See Eqs. (51) and (52)).
If the ∆-dependent cross section could be measured at a future electron ion collider for different values of x, and a similar systematic change of the dip position observed, it would be an interesting indication that we are approaching the saturated regime. Note that of course there are some caveats, since the detailed shape of the proton is not known and here we do not consider the growth of the proton with decreasing x [47][48][49][50], which will likely affect the detailed quantitative result.

VII. CONCLUSIONS
We have studied analytically and semi-analytically the properties of coherent diffractive dijet production in electron-proton and electron-nucleus collisions, using two different saturation models including impact parameter (and angular) dependence: the Golec-Biernat Wusthoff model and an Impact Parameter dependent McLerran Venugopalan model.
We derived general relations connecting angular correlations of the dipole orientation and impact parameter vector in coordinate space with angular correlations between mean dijet transverse momentum and hadron recoil momentum (Eq. (20)). We showed that the n−th Fourier harmonic of the amplitude for diffractive dijet production (in momentum space) depends only on the n−th harmonic of the dipole amplitude in coordinate space.
In the limit of large photon virtuality Q 2 and/or quark mass, the differential cross sections and elliptic anisotropies of the GBW model can be expressed completely analytically. In this limit, the P -and ∆dependencies factorize, providing distinct information on the projectile and target. The P -dependence showed interesting analytic structures such as a dip in the differential cross section (for longitudinal photon polarization), and changes in sign of the elliptic anisotropy (for both longitudinal and transverse polarizations), and provided insight into how these features depend on Q 2 , the quark mass, and the longitudinal momentum fractions of the quark and anti-quark. The ∆-dependence directly probes features of the target, being sensitive to the Fourier transform of the transverse density distribution.
In the case of the more realistic IP-MV model, where the anisotropy is explicitly driven by the gradients of the target geometry, we found approximate analytic expressions for the differential cross sections and elliptic anisotropies. In particular, we observed that the Pdependence was modified from the GBW model because of the presence of the logarithm in the dipole amplitude. Both the locations of dips in the longitudinal cross section and sign-change in the elliptic anisotropies shifted to larger values of P .
Approaching the saturation limit, we discussed the ex-pected modification to the features mentioned above; in particular their dependence on the saturation scale Q s . A more detailed analysis of the effects of approaching saturation could only be performed by numerical evaluation of our semi-analytic expressions. The numerical results confirmed several expectations. We observed an increase of the value of P at the dip position in the longitudinal differential cross section with increasing photon virtuality Q 2 and saturation scale Q s (increasing mass number A). A similar behavior was found for the change in sign in the elliptic anisotropy (for longitudinal photons). At low Q 2 the transversely polarized differential cross section and elliptic anisotropy dominate over their longitudinally polarized counterparts. The difference between them decreases with increasing Q 2 .
As a function of momentum transfer ∆, we observed the anticipated diffractive dips in the differential cross section for gold. We also observed a dip in the ∆dependent cross section for a proton target, despite the Gaussian shape of the assumed proton density profile. In the latter case, dips appear because of the unitarization of the dipole amplitude, signaling the effects of saturation, which leads to an effectively non-Gaussian shape of the proton. The ∆-dependence of the elliptic anisotropy also showed some effects of saturation such as the change in sign for the case of gold.
To gain more insight into the effects of saturation, we studied the x-evolution of the differential cross section at fixed values of P and ∆. We observed the expected slowdown in the growth of the differential cross section with decreasing x, with the effect setting in at larger values of x for the larger targets because of their larger saturation scales at a given x.
Finally, we studied the ∆-dependence of the diffractive dijet cross section for a proton target and different values of the saturation scale (representing varying x values), and observed a decrease of the value of ∆ at the diffractive dip position with increasing Q s . We compared to the "non-saturation" model, where the dipole was expanded to quadratic order, and did not observe any dips as expected. We argued that if such diffractive dips and their dependence on x could be measured experimentally in diffractive dijet events in e+p collisions, this could provide a strong indication of the presence of saturation effects.
The semi-analytic approach presented in this paper does not include some potentially important physical effects of the small x evolution, such as the growth of the color charge profile with decreasing x, and the corresponding modification of the color charge density gradients. In this work the x-evolution was only incorporated via the parametrization of the saturation scale Q s in Eq. (60). The effect was included in numerical studies using g Jalilian Marian-Iancu-McLerran-Weigert-Leonidov-Kovner evolution [51][52][53][54][55][56][57], and has been analyzed in [11]. The evolution of the dipole could also be studied using the BK evolution [58,59] with impact parameter dependence as studied in [60]. Also, a more detailed analysis will require the incorporation of the dependence of x on the dijet momenta [11].
Nevertheless, what makes our approach a very powerful tool for understanding what physical features of projectile and target are important for the process of diffractive dijet production in e+p and e+A collisions, is that we were able to find fully analytic expressions for cross sections and elliptic anisotropies in certain limits. Even the semi-analytic expressions only involve simple integrals that are easily evaluated numerically, which is especially helpful for examining the regimes of large P and/or ∆. This allows us to efficiently constrain the most interesting setup and kinematic regions in future experiments. In particular, we provided predictions for values of P and ∆ that maximize the magnitude of the elliptic anisotropies for different targets to assist future experiments in observing these interesting correlations.
Finally, we point out that the presented techniques could be extended to study other processes such as inclusive dijet production. In this case, one needs to work out an expression for the quadrupole, incorporating the effects of the geometry of the target. Similarly, one could attempt to extend this analysis to inclusive diffractive dijet and incoherent diffractive dijet production. In the latter case we expect to gain sensitivity to the local structures of the target and possibly angle dependent fluctuations. space to momentum space In this appendix we prove the relation in Eq. (20). We first lay out the conventions for the Fourier transforms and mode expansion. The Fourier transform and inverse Fourier transform are normalized as follows For rotationally symmetric functions, the Fourier mode decomposition is given by where the F k andF k can be computed by projection To prove Eq. (20), we use the following identity Then the Fourier transform Eq. (A1) can be expressed as a Bessel expansion. Changing the variables to θ rb and Θ = 1 2 (θ r + θ b ), we arrive at The Θ integral is trivial and is proportional to a Kronecker-Delta, 2πδ n,−m , which we use to contract the summation in m, yielding where we used J −n (z) = (−1) n J n (z).
To perform the angular integral, we plug in the Fourier mode expansion Eq. (A3) and find Contracting the summation in n and comparing with the expansion in Eq. (A4) we find Representation of Bessel functions of the first kind Representation of modified Bessel functions of the first kind The following integral is the backbone for our analytic computations: By taking derivatives with respect to P or ε f , and using recurrence relations for derivatives of J n (z) and K n (z), one finds It is also useful to have expressions for the Fourier transform of an isotropic function. They follow from the standard definition of the Fourier transform and Eq. (B1):T (∆) = 2π bdbJ 0 (∆b)T (b) , One can obtain interesting relations by taking derivatives. For example: where we used By inverting Eq.(B7), one has Appendix C: Details of analytic calculations of differential cross section and elliptic anisotropy In order to compute the differential cross sections and elliptic anisotropies (Eqs. (28) and (36)), it is enough to calculate the functions in Eqs. (26), (27), (34), and (35). For the sake of simplicity we ignore the small corrections to the differential cross section, i.e. we only keep the terms |F 0 | 2 and |∂ PG0 | 2 ). In this appendix we show the explicit calculations for these expressions in the limit Q s ε f , in which the dipole amplitude can be expanded to quadratic order. We start with the GBW model, for which we find exact analytic results, and then proceed to derive approximate expressions for the impact parameter dependent MV model.
The other expressions can be obtained in a similar fashion to read whereT 2 (∆) = 2π bdbJ 2 (∆b)T (b) is the 2nd order Hankel transform of T (b).

Impact parameter dependent McLerran Venugopalan model
We now consider the Impact Parameter dependent McLerran Venugopalan model in the limit Q s ε f . The expressions forF 2 and ∂ PG2 can be solved exactly. For example, one has F 2 (P, ∆) = 2π 8 Q 2 s rdrJ 2 (P r)r 2 K 0 (ε f r) Using Eqs. (B5) and Eq. (B8) to solve the r and b integrals, respectively, we obtaiñ Similarly, one has ∂ PG2 (P, ∆) = −2πQ 2 s P (P 2 − ε 2 f ) The expressions forF 0 and ∂ PG0 , on the other hand, cannot be solved exactly due to the presence of the logarithm in the r dependent part of the integrand. For example, one has F 0 (P, ∆) = π 2 Q 2 sT (∆) rdrJ 0 (P r)f 0 (r) , with f 0 (r) = r 2 log 1 m 2 r 2 + e K 0 (ε f r) . (C10) It would be useful to approximate Eq. (C10) by an expression of the form of r 2 K 0 (ε f r) as it appears in the GBW model, for which we had an analytic solution. First, one should note that the convolution (Fourier transform) in Eq. (C9) is dominated by the maximum of f 0 (r). Thus, in the following we focus on reproducing the effect of the modified location and height of the maximum of Eq. (C10). Because of the logarithmic factor, Eq. (C10) develops a maximum at a smaller value of r compared to r 2 K 0 (ε f r), which depends on the ratio κ ≡ ε f /m. We will assume that κ 1. To see this more explicitly, we change to the variable u = ε f r f 0 (u) = u 2 ε 2 f log κ 2 u 2 + e K 0 (u) .
The maximum of this function occurs at (ignoring the factor of e inside the logarithm) while in the GBW model the maximum occurs at u max = 2K 0 (u max )/K 1 (u max ) ≈ 1.5. Therefore, we see that in the IP-MV model, the location of the maximum is shifted to a smaller value of u (compared to GBW): with ξ = log(κ 2 ) log(κ 2 )−1 . For values of κ = 3 − 10, one has ξ= 1.3 -1.8.
We thus approximate f 0 in Eq.(C10) by f 0 (r) ≈ C 1 r 2 K 0 (ξε f r) , where C 1 = log κ 2 K 0 (1.5/ξ)/K 0 (1.5) > 1. This expression reflects the shift in the location of the maximum and the increase in the height of the maximum.
transfer kicks q ∼ 1/R with each scattering. The approximation above then is valid if 1/R < 2m. For a large target such as a nucleus this is satisfied, while for a proton the approximation is questionable. Since the single scattering momentum transfers are restricted to q 2m, we will not trust this approximation much beyond ∆ ∼ 2m. Even though one might be concerned about the divergence of as m → 0 in Eq. (43), one should note that in the lim m→0 Θ(q, m) = 1/2; whose effect is to replace the regulator m by the finite system size, or 1/R in Eq. (43).

Appendix E: Nuclear saturation scale
The local saturation scale Q 2 s (b) = Q 2 s T (b) is proportional to the charge density squared of the target at point b. Since the nucleons are assumed to be uncorrelated, the total charge squared in the nucleus is the sum of the charges squared of all its nucleons. Thus, one has For our choice of proton profile Eq. (57) and nuclear profile Eq. (58), and using where we assumed R A a A to approximate the ρ A in Eq. (59) by a hard sphere and N A ≈ 1/(2R A ). Thus we have where we used R p = 0.4 fm and the approximate expression R A = 1.1A 1/3 fm for large nuclei. A similar expression was obtained in [61], where the authors assumed a cylindrical shape for nuclei. An expression assuming spherical nuclei and nucleons was obtained in [62].
If one accounts for the non-zero a A , then one finds the following relation between the saturation scales: Q 2 sAu = 2.17 Q 2 sp and Q 2 sCu = 1.44 Q 2 sp , for gold and copper respectively.