Diffractive Dijet Production and Wigner Distributions from the Color Glass Condensate

Experimental processes that are sensitive to parton Wigner distributions provide a powerful tool to advance our understanding of proton structure. In this work, we compute gluon Wigner and Husimi distributions of protons within the Color Glass Condensate framework, which includes a spatially dependent McLerran-Venugopalan initial configuration and the explicit numerical solution of the JIMWLK equations. We determine the leading anisotropy of the Wigner and Husimi distributions as a function of the angle between impact parameter and transverse momentum. We study experimental signatures of these angular correlations at a proposed Electron Ion Collider by computing coherent diffractive dijet production cross sections in e+p collisions within the same framework. Specifically, we predict the elliptic modulation of the cross section as a function of the relative angle between nucleon recoil and dijet transverse momentum for a wide kinematical range. We further predict its dependence on collision energy, which is dominated by the growth of the proton with decreasing $x$.


I. INTRODUCTION
Diffractive processes in deep inelastic scattering of electrons off protons or heavier nuclei can provide important information on the target's structure in coordinate and momentum space [1]. For example, diffractive vectormeson production is sensitive to the spatial profile and fluctuations in the target for the coherent and incoherent cross sections, respectively [2][3][4][5][6].
A future electron-ion collider (EIC) [40][41][42] will allow proton tomography with unprecedented precision, measuring parton position, momentum and spin inside protons and nucleons. The EIC would open up vast possibilities for understanding the gluon Wigner distribution by measuring diffractive dijet production as a function of the nucleon recoiled momentum and the dijet transverse momentum. Similar research directions may be explored at a possible Large Hadron Electron Collider [43].
At very small x the color glass condensate effective theory, describing quantum chromodynamics (QCD) in the high-energy limit, is a suitable framework to compute diffractive processes. In this framework, slow modes in the fast-moving target are highly occupied gluon fields, which can be described classically by the Yang-Mills equations. Fast modes act as sources for these slow modes. Renormalization group equations govern the evolution of the separation between hard and soft modes towards lower momentum fractions. These are the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equations [82][83][84][85][86][87][88], which in the limit of a large number of colors simplify to the Balitsky-Kovchegov equations [79,80]. The JIMWLK equations resum all terms enhanced by large logarithms in 1=x, resulting in leading logarithmic accuracy. Next-to-leading-order corrections were included several years ago [89][90][91][92], but, as opposed to the leading logarithmic case [93,94], at this point no (numerical) solution of the NLO JIMWLK equations has been explored.
The CGC provides the necessary ingredients to compute coherent diffraction in the dipole picture. Here, colorsinglet Pomeron exchange is understood as the colordiagonal [2] interaction of a dipole with a stochastic ensemble of classical color fields in the eikonal limit. A stochastic average over target color configurations is required, which when performed on the level of the scattering amplitude is equivalent to assuming that the target remains intact ("coherent diffraction").
In the dipole picture, hard diffractive dijet production in small-x DIS has been considered in Ref. [7] and in Ref. [95]. Related work includes Ref. [96], where it was suggested to study forward diffractive quarkonia production in p þ p collisions to probe the Weizsäcker-Williams gluon distribution and Ref. [97], where access to the gluon Wigner distributions via ultraperipheral p þ A collisions was discussed. Exclusive double production of pseudoscalar quarkonia and its relation to the gluon Wigner distributions in nucleon-nucleon collisions was explored in Ref. [98]. Exclusive diffractive two-and three-jet production in photon-hadron scattering within k T factorization and the cancellation of IR, collinear and rapidity singularities for the two-jet cross section was studied in Ref. [99].
In this manuscript, we present calculations of the gluon Wigner and Husimi distributions in the proton within the CGC framework at leading logarithmic order. We introduce a spatially dependent color charge distribution of the proton, constrained by DIS data from HERA. The energy dependence of the corresponding Wigner distribution is determined by numerically solving the JIMWLK evolution equations. We compute correlations between the impact parameter and transverse momentum and extract the elliptic anisotropy coefficient of Wigner and Husimi distributions.
In order to connect the Wigner distribution to experimental observables, we perform an extensive study of diffractive dijet production cross sections in e þ p collisions at typical EIC kinematics within the same framework. We focus on dijet kinematics in the so-called correlation limit and predict the dependence of elliptic modulations as a function of the relative angle between the nucleon recoil and dijet transverse momentum on collision energy and kinematics, which can be tested at a future EIC.
This manuscript is organized as follows. In Sec. II A we illustrate how gluon Wigner distributions can be directly computed within the CGC in the small-x limit. In Sec. II B we discuss the production cross section for dijets in virtual photon-nucleus scattering, for arbitrary photon virtuality and quark mass. In Sec. III we compute gluon Wigner and Husimi distributions from the CGC and study correlations between the impact parameter and transverse momentum.
In Sec. IV we present our results for diffractive dijet production cross sections in e þ p collisions at typical EIC energies. First, in Sec. IVA we present a simple baseline study based on the impact parameter-dependent saturation (IP-Sat) model and disentangle genuine correlations from kinematic effects. We then compute diffractive cross sections for charm jets from the CGC in Sec. IV B and study the dependence of the elliptic Fourier coefficient on the dijet momentum, photon virtuality and collision energy.
Our findings are supplemented by multiple appendices. In Appendix A we present details of the Wigner function derivation from the CGC and a practical approach to numerically compute Wigner and Husimi distributions in this framework. In Appendix B we discuss conventions for transverse dijet momentum variables and in Appendix C we investigate the dependence of azimuthal dijet correlations on the proton size.
Most inclusive processes are not very sensitive to the full momentum and spatial structure of the nucleon, as the impact parameter and/or transverse momenta are usually integrated over. In contrast, more exclusive cross sections can provide access to five-dimensional gluon Wigner distributions and diffractive dijet production is particularly well suited for this task [7]. To probe impact-parameterand transverse-momentum-dependent gluon dipole distributions, the authors of Ref. [7] proposed to study coherent diffractive dijet production off a nuclear target in the correlation limit (defined below), a process which can be probed at the future electron-ion collider.

B. Diffractive dijet production in the dipole picture
In this section, we outline the basic elements for the calculation of coherent diffractive dijet production in DIS off a nuclear target, where l, l 0 (P, P 0 ) are in-and outgoing lepton (target) fourmomenta, q ≡ l 0 − l and p 0 ≡ ðp þ 0 ; p − 0 ; p 0 Þ and p 1 ≡ ðp þ 1 ; p − 1 ; p 1 Þ are the four-momenta of the produced dijet. In particular, p þ 0 ≡ zq þ , p þ 1 ≡zq þ denote the light-cone momenta of the jets, where q þ is the photon longitudinal momentum andz ¼ 1 − z.
The scattering process occurs via the exchange of a virtual photon between the target and a lepton, and can be interpreted as photon-nucleon scattering for given lepton kinematics, γ Ã ðqÞ þ NðPÞ → N 0 ðP 0 Þ þqðp 0 Þ þ qðp 1 Þ. In the dipole approximation, the scattering then entails the virtual photon fluctuating into a quark-antiquark color dipole, which in turn interacts with the color field of the nucleus. In the eikonal approximation the quarks exchange color with the target and receive a "kick" of transverse momentum. In coherent diffractive scattering no net color is exchanged between the projectile and target and the target stays intact. Below, we call outgoing quarks and antiquarks "jets," but a fully phenomenological study requires fragmentation via event generators [132].
Following Refs. [95,133], we write the S-matrix element for the diffractive process γ Ã N →qqN 0 as where q f ¼ Z f e, f 0;1 , h 0;1 ; a; b are the charge, flavor, helicity and color of the outgoing dijets, UðxÞ are the Wilson lines in the fundamental representation (discussed below) and Here, Ψ λ ðr; Q 2 Þ is the photon light-cone wave function with longitudinal or transverse polarization λ ¼ 0; AE1 describing the splitting into a dipole of size and orientation r ¼ x 0 − x 1 [133,134].
For the coherent diffractive process, the target remains intact and we average over target color configurations on the amplitude level. The production cross section of a dijet pair with momenta (p þ 0=1 , p 0=1 ) can be written as where M is the scattering amplitude, defined as where S is the S-matrix. A diffractive process is characterized by color neutral exchange ("Pomeron exchange") and the resulting rapidity gap between the outgoing dijet system and the struck target, where no particles are produced. The scattering matrix is color diagonal and allows us to write the dipole-target interaction as [2] In this study, we express dijet production cross sections in the transverse momentum variables 3 In these coordinates, the dijet production cross section for transversely polarized photons reads [95,133] where N is given by Eq. (4) and we abbreviated Q ≡ ðzzQ 2 þ m 2 q Þ 1=2 , α EM ¼ e 2 =4π and ζ 2 ≡ z 2 þz 2 . Here, ϵ λ ðqÞ ¼ 1= ffiffi ffi 2 p ð1; λiÞ are light-cone photon polarization vectors and K 0;1 are modified Bessel functions of the second kind. The longitudinal momentum fraction of the "Pomeron" is defined by where Q 2 ≡ −q 2 ¼ −ðl − l 0 Þ 2 is the virtuality of the photon, W 2 ¼ ðP þ qÞ 2 is the center-of-mass energy squared of the photon-target system, m N is the target mass, M 2 is the invariant mass-squared of the dijet system and t ≡ −ðP 0 − PÞ 2 ¼ −Δ 2 is the Mandelstam variable. Analogously, the cross section for photons with longitudinal polarization reads Angular correlations can be parametrized by azimuthal Fourier decomposition, Here, θðP; ΔÞ ≡ θðPÞ − θðΔÞ ð 14Þ is the relative angle between the dijet and target recoil momentum. The coefficients v n are experimentally measurable by analysis of the azimuthal distribution of dijets. By Fourier transformation, the angular correlation in the dijet cross section is sensitive to the relative orientation between r and b.
Here, we study the dijet cross section in the correlation limit, jPj ≫ jΔj, where the individual jets are almost back to back. In this limit, one is most sensitive to dipole contributions from the peripheral region at large jbj. In Ref. [97] it was further shown that in the correlation limit for Q 2 ¼ 0, a direct relation between the diffractive dijet cross section and the gluon Wigner distribution can be established.
C. The dipole amplitude at small x P In this section, we discuss two approaches to compute the dipole amplitude, Eq. (4). The first will be a model parametrization from the IP-Sat model [108], which has been successfully used to describe a range of data, from HERA inclusive and diffractive e þ p DIS data [110,112,135] to n-particle multiplicity distributions in p þ p and p þ A collisions at RHIC and LHC [113,136,137]. The second one will be a CGC computation, where the dipole interaction with the target is directly encoded in fundamental Wilson lines.

IP-Sat model
In the IP-Sat model the dipole amplitude is given by Its impact parameter dependence arises from the transverse spatial color profile of the proton, assumed to be Gaussian, The IP-Sat model is x P -dependent through the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [138][139][140][141] evolved gluon distribution function x P gðx P ; μ 2 Þ at the scale The proton width B p , the initial scale μ 0 and the initial conditions for x P gðx P ; μ 2 Þ are obtained from fits to HERA DIS data [112]. Note that, independent of jbj, the IP-Sat model always parametrizes the long-distance behavior as N → 1 for jrj → ∞, in contrast to the CGC computation for finite systems [107,142].
The IP-Sat model is not useful when studying azimuthal correlations in diffractive dijet production, because it does not depend on r · b, but only on r 2 and b 2 separately. Consequently, it cannot reproduce any angular correlation between P and Δ. In addition, the IP-Sat model contains only DGLAP evolution and we expect it cannot capture important features of the full JIMWLK CGC computation. However, by artificially introducing angular correlations, we can consider the IP-Sat parametrization as a simple baseline test, where we included a θðr; bÞ ≡ θðrÞ − θðbÞ-dependent term Varying the parameterc allows one to regulate the amount of anisotropy in the gluon Wigner distribution by hand and to study the effects on dijet cross sections. A similar parametrization was used in Ref. [95], and a more sophisticated analytic expression, albeit with no x P dependence, was discussed for example in Ref. [143]. In the latter study, in accordance with our expectations, cos½2θðr; bÞ correlations occurred in connection with gradients of the target color charge density. Specifically, these correlations vanish in the homogeneous regime at small b and r, a feature not reproduced by our Eq. (17). Instead of performing further modifications of the IP-Sat model, we set out to compute the gluon dipole distribution directly from the color glass condensate effective theory at small x P , including its energy evolution via the JIMWLK renormalization group equations in the next subsection.

Color glass condensate computation
In this section, we discuss the derivation of the dipole amplitude from the color glass condensate effective theory [74][75][76][77][78]. For an initial x P , it is computed as a stochastic average from fundamental Wilson lines UðxÞ, which are obtained by solving classical Yang-Mills equations with target color sources sampled from a local Gaussian distribution where the color charge density is related to the IP-Sat value for the saturation scale at moderately large x P , An impact-parameter-independent target allows to compute the value of the parameter c directly [144], and we will vary the parameter between 0.75-0.85 in the impact-parameterdependent case. For a given target color configuration, the solutions of the Yang-Mills equation specify Wilson lines, Here, an IR cutoffm ∼ 0.2-0.4 GeV has been introduced to avoid nonphysical Coulomb tails at distances ∼Λ −1 QCD where nonperturbative effects become important.
After computing Wilson lines at some initial x P , we determine the energy evolution of the dipole amplitude using the JIMWLK equations [82][83][84][85][86][87][88]. At leading logarithmic order, the JIMWLK renormalization group equations can be written as a functional Fokker-Planck equation [93] and explicitly expressed as a Langevin equation for the stochastic Wilson lines UðxÞ, where dy ¼ dx P =x P is the rapidity and t a are SUðN c Þ generators in the fundamental representation. The drift term in Eq. (22) reads where the Wilson linesŨðzÞ,Ũ † ðxÞ and generators T a are in the adjoint representation. The ξ b i ðz; yÞ's are stochastic and are sampled from a Gaussian distribution with zero mean and variance given by The kernel in Eq. (22) is where we abbreviated It is possible to eliminate the drift term by writing [145] Within the CGC JIMWLK evolution, long-distance tails are encountered in the Coulomb kernel K i . These lead to an exponential growth of the cross section with rapidity [107,146,147], ultimately violating the Froissart unitarity bound [148,149] unless regulated by nonperturbative physics at large distance scales. To regularize the nonperturbative regime, we follow the prescription of Ref. [142] and instead use the kernel where m ∼ Λ QCD and K 1 is a modified Bessel function. We include running coupling effects in our analysis and evaluate the coupling constant as follows [145]: where r ≡ jx − zj, μ 0 ¼ 0.28 GeV and χ ¼ 0.2 [6]. In practice, we fix the initial condition for the JIMWLK evolution at x P ¼ 0.01 from the IP-Sat model. We then evolve towards smaller x P by solving the JIMWLK evolution equations. We study different choices of IR regulatorsm and m, in a range which is constrained by and consistent with data [5,6], and we compare fixed vs running coupling α s , cf. Eq. (28). For more details, we refer the reader to Ref. [6].

III. GLUON WIGNER AND HUSIMI DISTRIBUTIONS AT SMALL x FROM THE COLOR GLASS CONDENSATE
In the CGC approach, angular correlations between the impact parameter and dipole orientation in the dipole amplitude N ðr; b; xÞ emerge naturally and need not be modeled, in contrast to e.g., the situation in the IP-Sat model (17).
In this section, we present a CGC study with fixed α s ¼ 0.21 and infrared regulatorsm ¼ 0.4 GeV and m ¼ 0.2 GeV as in Ref. [6], and obtain the dipole amplitude by solving the JIMWLK evolution equations with the rapidity defined as where xð0Þ ¼ 10 −2 . In Fig. 1, we show the normalized dipole amplitude as a function of the relative angle θðr; bÞ between the impact parameter b and dipole orientation r.
We present results for the initial condition (y ¼ 0) and after y ¼ 1.5 and y ¼ 3 units of rapidity evolution. Here, the dipole amplitude is largest whenever the dipole is oriented along the impact parameter, which is expected as this configuration is most sensitive to (radial) color gradients (for a similar analysis using impact-parameter-dependent BK evolution, see Ref. [147]).
To quantify this behavior, the elliptic component v 2 of the dipole amplitude is defined as where v 0 is the average dipole amplitude. We plot the x dependence of v 2 for different impact parameters in Fig. 2 and find that the energy (rapidity or x) evolution suppresses the elliptic component significantly. This is a consequence of the proton's growth with energy, leading to smoother density gradients for fixed impact parameter. We observe a similar trend for the energy dependence of diffractive dijet cross sections presented in Sec. IV B. Further discussion of the proton size dependence can be found in Appendix C. Next, we study the Wigner and Husimi gluon distributions, which encode information about the transverse momentum and coordinate dependence of the small-x gluons, as discussed in Sec. II A. The quasiprobabilistic gluon Wigner distribution xWðx; P; bÞ is real but not necessarily positive; following Ref. [23] its small-x (CGC) limit can be computed from Eq. (3), and reads The Husimi distribution is positive semidefinite and can be computed from the Wigner distribution by Gaussian smearing of the transverse momentum and impact parameter. Following Ref. [23], it can be written as where l is an arbitrary smearing parameter. We investigated various values for l and present results for l ¼ 1 GeV −1 only. This is a reasonable choice, as it ensures that the spatial smearing scale is smaller than the proton, but also keeps the (inversely proportional) momentum smearing scale small enough to resolve the transverse momentum spectrum of soft gluons in the proton's wave function. However, the dependence on this smearing scale is a disadvantage for interpretations of the Husimi distribution. We argue below that while Husimi and Wigner distributions agree in certain limits where smearing has no effect, some of the former's features can be l dependent, making it more favorable to work with the Wigner distribution instead.
Similarly as for the dipole amplitude, we study the elliptic modulation of Wigner and Husimi distributions.
Because v 2 [in Eq. (30)] is not defined at zero crossings of the respective distribution (recall that the Wigner distribution is not positive definite), we proceed with a different parametrization [23] and similarly for the Husimi distribution, xHðP; b; xÞ ¼ xH 0 þ 2xH 2 cosð2θðP; bÞÞ: ð34Þ In Fig. 3, we show the lowest moments xW 0 and xH 0 of Wigner and Husimi distributions as a function of transverse momentum and for different rapidities y ¼ 0 and y ¼ 1.5. The effect of smearing is most prominent at small momenta, where jPjl ≪ 1. At large momenta, where both Husimi and Wigner distributions are positive, smearing has no effect and the distributions agree. The presented Wigner distribution should be contrasted with the results of a previous analysis of gluon Wigner distributions [23], with results relatively similar to ours. Overall we find qualitative agreement for jPj ⪆ 0.1 GeV; however we do not find that the position of the characteristic peak in Fig. 3 evolves with energy in our study. We also do not recover a feature, seen in Ref. [23], where the Wigner distribution rises from below to zero at very small momentum.
We are not surprised by the differences between the model study [23] and our numerical CGC computation. At small momentum, (nonperturbatively) large dipoles (as large as a few fm) give a significant contribution to the  These results are averaged over the azimuthal angle between the transverse momentum and the impact parameter. 4 In Appendix A we outline a practical computation of these coefficients.
Wigner distribution, while the only effective regulator of these contributions is the proton size. In our framework this extreme IR behavior is parametrized by the infrared regulators m andm [see Eqs. (21) and (27)], which are constrained to some degree by HERA data. 5 In contrast, the authors of Ref. [23] introduced an explicit exponential suppression factor (by hand) to cut off contributions from this regime. 6 In Fig. 4 we show the elliptic components xW 2 and xH 2 as a function of transverse momentum. We find the magnitude of xW 2 to be on the percent level when compared to xW 0 , similar as in Ref. [23]. We observe that the elliptic part of the Wigner distribution is much larger than that of the Husimi distribution for jPj ≲ 2 GeV. It is numerically difficult to extract xH 2 at larger momentum, but we find that the distributions agree in this regime within error bands shown in Fig. 4. The small-jPj behavior of xW 2 and xH 2 is again sensitive to infrared regularization and deserves further study.
In order to gain intuition and to illustrate some of the difficulties in its interpretation, we study the Husimi distribution in more detail. For better comparison with the v 2 component of the dipole amplitude and of the dijet cross section studied below, we present the normalized elliptic coefficient in Fig. 5, Here, we show the momentum dependence of v H 2 for y ¼ 0 (x ¼ 10 −2 ) at different impact parameters. As expected, v H 2 vanishes at small jbj where the average proton color profile is nearly homogeneous and cannot have any angular dependence. The elliptic component peaks at a characteristic moderately large jPj (the exact position of the peak depends on the parameter l) and approaches zero at small and large momentum.
In Fig. 6, we show the energy (rapidity) dependence of v H 2 at fixed impact parameter, where again the growth of the proton with energy leads to a decrease of v H 2 , which may be compared to the v 2 component of the dipole amplitude and of the dijet cross section studied below. Remarkably, the decrease of jv H 2 j with energy, caused by the decreasing density gradients, is not uniform for all jPj. In the smallmomentum regime, the decrease of jv H 2 j is prevented by the fact that when the proton grows, one first starts to include dipoles with transverse separation jrj ∼ jPj −1 with large elliptic modulation. Eventually, the proton grows larger than these dipoles and the probed density gradients become smaller, resulting in a (delayed) decrease of jv 2 j at small transverse momentum.
This behavior exemplifies some of the difficulties in assigning a physical interpretation to the Husimi distribution. For example, with smaller l (less smearing in transverse coordinate space), the dipole sizes would be more strongly limited and jv H 2 j would uniformly FIG. 6. Rapidity evolution of the elliptic part of the Husimi distribution as a function of transverse momentum jPj. 5 A more detailed discussion of the IR regulator dependence can be found in Sec. IV. 6 Note that large distance contributions to the Husimi distribution are suppressed by the smearing scale l. decrease for all jPj. 7 While the positivity of the distribution might be an appealing argument for a probabilistic interpretation, we wish to argue that more robust information can be extracted from the Wigner distribution.
In the next section we compute diffractive coherent dijet production cross sections in electron-proton scattering at typical EIC energies. Here, we study angular correlations between the dijet transverse momentum and proton recoil, which are directly sensitive to correlations between the impact parameter and transverse momentum in the gluon Wigner and Husimi distributions.

IV. COHERENT DIFFRACTIVE DIJET PRODUCTION
In this section, we present results for coherent diffractive dijet production in virtual photon-proton scattering. We investigate angular correlations between transverse dijet momenta and target recoil. We focus on typical EIC kinematics and for most of our study we fix the proton beam energy E p ¼ 250 GeV, and the center-of-mass energy to be W ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðP þ qÞ 2 p ¼ 100 GeV. Our results are presented in the analysis frame, where the photon and incoming nucleon have no transverse momentum and our convention is that the photon (proton) has a large þ (−) momentum. Below, we first discuss results from the IP-Sat parametrization, which we use to disentangle kinematic effects from genuine correlations in the underlying gluon dipole distribution.

A. Baseline study: Angular correlation in the modified IP-Sat model
We begin with an overview of our results for coherent diffractive dijet production from the IP-Sat parametrization, with [Eq. (17)] and without [Eq. (15)] angular correlations between the impact parameter and dipole orientation. We present results for typical EIC kinematics, where E p ¼ 250 GeV, W ¼ 100 GeV and Q 2 ¼ 1 GeV 2 .
In Fig. 7 we show the cross section for longitudinally and transversely polarized photons as a function of transverse dijet momentum jPj, where we have set jΔj ¼ 0.1 GeV and θðΔ; PÞ ¼ π. These results are for light flavor jets with quark mass m q ¼ 0.03 GeV and symmetric longitudinal jet momenta z ¼z ¼ 0.5. As expected, the cross section is steeply falling with increasing jet momentum jPj, and the total cross section is dominated by the transversely polarized photons at small jPj and by longitudinal photons at larger jPj.
In Fig. 8, we show normalized dijet cross sections as a function of the relative angle θðΔ; PÞ between the jet transverse momentum P and target recoil Δ. Here, we integrate over the longitudinal momentum fraction z;z ∈ ½0.1; 0.9. We compare results from the conventional IP-Sat parametrization, Eq. (15), without angular correlations between the impact parameter and dipole orientation (thin lines), with those obtained from the modified IP-Sat model (17), including angular correlations (thick lines). Here, we have setc ¼ 1, which induces a nonzero dijet v 2 of the order of a few percent. The elliptic correlation has opposite sign for longitudinal (v 2;L < 0) and transverse photons (v 2;T > 0).
In Fig. 9 we present the longitudinal, transverse and total v 2 as a function of the dijet momentum jPj with kinematics as in Fig. 8 The cross section is integrated over z ∈ ½0.1; 0.9 and over the direction of the proton recoil θðΔÞ. 7 One could optimize the choice of l for each y, thereby finding the optimal resolution for a given proton size and characteristic transverse momentum scale Q s . In this way one could reduce the parameter dependence of the elliptic modulation in the Husimi distribution, but this is not a very attractive option. transverse components have opposite sign and the total v 2 is partially cancelled between the two contributions. Figure 10 shows the normalized transverse (left), longitudinal (center) and total (right) dijet cross sections as a function of θðΔ; PÞ, separately for different values of z. Remarkably, a nonzero v 1 ðzÞ ≡ hcos½θðΔ; PÞi z is observed (most visible in the left panel) whose origin we discuss below. To illustrate this feature, we plot v 1 ðzÞ and v 2 ðzÞ in Figs. 11 and 12 for transversely (dashed) and longitudinally polarized (dash-dotted) photons with (thick lines) and without (thin lines) correlations between the impact parameter and dipole orientation. The nonzero v 1 ðzÞ is independent of the angular correlation in Eq. (17). Its emergence is due to the energy dependence of the dipole amplitude [Eq. (11)]. To show this, we Taylor expand the dipole amplitude, close to but away from the asymptotic correlation limit to linear order in jΔj=jPj, The linear term can be computed from Eq. (11), If we insert Eq. (37) into Eq. (36) and use it to compute the dijet cross section, the latter must have a nonzero v 1 ðzÞ at OðjΔj=jPjÞ, vanishing only in the exact correlation limit. However, because v 1 ðzÞ ¼ −v 1 ðzÞ, v 1 ðzÞ can be eliminated by z-symmetric integration; see for example FIG. 10. Cross sections for transversely (left) and longitudinally (center) polarized photons and the total cross section (right), as a function of the longitudinal momentum fraction of one jet z ¼ p þ 0 =q þ . A sizable v 1 is induced by kinematic effects, because x P is not constant along θðΔ; PÞ for fixed z, according to Eq. (11). This effect is OðjΔj=jPjÞ and vanishes only in the asymptotic correlation limit.   Figure 12 shows that a finite v 2 only appears forc ≠ 0. For longitudinal polarization v 2 is negative for all z and jv 2 j is maximal at z ¼ 0.5. In contrast, for transverse polarization, v 2 is positive for all z and maximal in the most asymmetric cases z → 0 and z → 1.
Next we compute the dijet cross section for the case of charm quarks. A simple analysis of Eq. (10) and Eq. (12) shows that dipoles with size jrj ⪆ Λ −1 QCD contribute to the dijet cross section for light quarks and small photon virtuality of only a few GeV 2 . While the IP-Sat dipole parametrization reports this limit as N ðjrj → ∞Þ → 1, no reliable theoretical description of angular correlations exists in this regime. Here, nonperturbative hadronic effects must be taken into account. From now on we focus on charm dijets, which by means of the large mass scale, allows us to predict cross sections for photon virtuality down to Q 2 ∼ 0 without unconstrained contributions from large distances ∼Λ −1 QCD . Results for the charm dijet cross section from the modified IP-Sat model are shown in Fig. 13. Here, we plot the jPj dependence of the transverse, longitudinal and total diffractive charm-dijet production cross sections for charm quark mass m c ¼ 1.28 GeV, photon virtuality Q 2 ¼ 1 GeV 2 , jΔj ¼ 0.1 GeV, θðΔ; PÞ ¼ π and z ¼z ¼ 0.5. In the shown range of jPj the transverse charm dijet cross section is significantly larger than the longitudinal component. A novel feature is the sharp drop of the longitudinal cross section around jPj ∼ 1.5 GeV. A similar, but less drastic, dip is visible for the transverse cross section.
These minima reflect sensitivity to the size of the projectile jrj and are understood from Fourier transformation to momentum space. A simple estimate of the characteristic inverse "size" of the photon from its wave function in coordinate space yields jr γ j −1 ∼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m 2 c þ zzQ 2 p ≈ 1.4 GeV, which roughly coincides with the minima of Fig. 13. 9 In Fig. 14, we show the jPj dependence of the transverse and longitudinal v 2 for charm jets from the modified IP-Sat  8 We note that a different choice of dijet and target recoil transverse momenta,P ≡zp 0 − zp 1 andΔ ¼ p 0 þ p 1 , discussed in Appendix B, avoids the cos½θðΔ; PÞ dependence [132]. These variables cannot be interpreted as Fourier conjugates to the impact parameter and dipole orientation and other kinematic effects must be considered. 9 The interpretation of these features is thereby analogous to the origin of minima usually observed in the t ¼ −jΔj 2 spectrum of diffractive cross sections. Here, diffractive minima indicate the size of the target, not the projectile. We also note that the lightquark photon wave function is "larger" than that for charm quarks where we expect a similar feature below 1 GeV (not shown in Fig. 7). We emphasize that the photon is "nonperturbatively large" for light quarks and small virtuality. model withc ¼ 1. We keep the kinematics as in the light quark case, jΔj ¼ 0.1 GeV, Q 2 ¼ 1 GeV 2 , z ∈ ½0.1; 0.9. A consequence of the minima shown in Fig. 13 is a strong increase of the longitudinal v 2;L at jPj ∼ 1.5 GeV where it changes sign, with v 2;L > 0 for jPj ≲ 1.5 GeV and v 2;L < 0 for jPj ≳ 1.5 GeV. The transverse component, too, changes sign which we attribute to the relative importance of the two terms in Eq. (10). Specifically, in Eq. (10) the second, mass-dependent term dominates for large jPj. Because this term behaves similarly as the corresponding contribution to the longitudinal cross section, it also results in a negative overall v 2 .
In the following section, we compute the dipole amplitude directly from the color glass condensate effective theory where angular correlations are included ab initio and we will reliably extract the energy dependence of our results.

B. CGC computation
In this section, we compute dijet cross sections from the CGC, as outlined in Sec. II C 2. For this quantitative study, we do not show light quark results because of the aforementioned sensitivity to large jrj ∼ Λ −1 QCD contributions which are not under good control theoretically.
We plot transverse, longitudinal and total cross sections in Fig. 15, where jΔj ¼ 0.1 GeV, Q 2 ¼ 1 GeV 2 , z ¼z ¼ 0.5 and θðΔ; PÞ ¼ π. We vary the infrared regulators m;m [see Eqs. (21) and (27)] and the value of the constant coupling α s , as well as the parametrization with a running coupling (28). We match the value of the saturation scale according to Eq. (20) with c ¼ 0.75 for an IR regulator m ¼ 0.4 GeV, and c ¼ 0.85 form ¼ 0.2 GeV, thus ensuring a consistent overall dipole normalization. We further adjust the value of the fixed α s , when we change m in the JIMWLK kernel (27) to ensure comparable energy evolution speed. We use this specific parameter range, because it has been used to compute proton structure functions and diffractive vector meson cross sections in Ref. [6], and was shown to be consistent with HERA data. Shown in Fig. 15, the CGC-parameter dependence is very small for charm-dijet cross sections, indicating a negligible sensitivity to large distances jrj ∼ Λ −1 QCD . For comparison we include the IP-Sat results (gray dashed curve), which differ significantly from the CGC in magnitude, but show similar qualitative features. We attribute the different normalization to the fact that in the CGC the small-jrj regime of the dipole amplitude carries a larger relative weight, compared to IP-Sat, as shown in Fig. 5 of Ref. [6]. For charm dijets the relevance of this small-jrj regime is enhanced.
In Fig. 16 we plot the v 2 obtained from the CGC computation for transverse, longitudinal, and total cross sections as a function of jPj, integrated over azimuthal angle θðΔÞ and z;z ∈ ½0.1; 0.9. Figure 16(a) shows the transverse v 2;T for all CGC parameter sets discussed above. All parametrizations agree qualitatively and the results are robust under variation of the JIMWLK evolution parameters (stars, circles, and triangles). Including the running coupling α s via Eq. (28), vs a fixed coupling only minimally effects the JIMWLK evolution. Results with different IR regulatorsm differ by up to 20% (squares vs stars, circles, and triangles). Both valuesm ¼ 0.2 GeV and m ¼ 0.4 GeV are in agreement with HERA data [6] and further comparison with data is necessary to reduce this systematic uncertainty.
In Fig. 16(b) we show results for the longitudinal anisotropy coefficient v 2;L , where the most prominent feature is the rapid growth of v 2;L in the jPj region where  , and total (c) cross sections for charm dijets from the CGC for the longitudinal cross section drops; see Fig. 15(b). All CGC parametrizations give qualitatively and quantitatively similar results, and we attribute the largest uncertainty to the dependence onm. The total v 2 , which is a weighted sum of transverse and longitudinal components, is shown in Fig. 16(c). At small jPj, longitudinal and transverse cross sections have positive v 2 and the total cross section is dominated by the transverse component. At large jPj ∼ 2.5 GeV, the longitudinal v 2;L is negative while again the total v 2 is dominated by photons with transverse polarization.
In Fig. 17 we study the energy dependence of v 2 , by varying the center-of-mass energy W, while keeping the dijet kinematics fixed (at jPj ¼ 1 GeV, Δ ¼ 0.1 GeV, with photon virtuality Q 2 ¼ 1). We integrate over z ∈ ½0.1; 0.9 and over the azimuthal angle θðΔÞ. Shown are results for three different energies W ¼ 60, 100, 150 GeV. Because our results are integrated over z and θðΔ; PÞ, we plot them as a function of the average x P in the respective integration range z ∈ ½0.1; 0.9 and θðΔ; PÞ ∈ ½0; 2π [see Eq. (11)], where horizontal bars denote the standard deviation of x P in this range. We compare results from the full CGC computation form ¼ 0.4, m ¼ 0.2 with fixed α s ¼ 0.21 (squares) to results without JIMWLK evolution (circles), where the x P dependence is entirely determined from the IP-Sat parametrization of HERA data. In the latter case, only the saturation scale Q s changes with energy, while in the full CGC computation the transverse size of the proton increases as well. This growth is accompanied by a relative change in the gradients of the target color density. Consequently, we expect a decrease of v 2 with energy, as is in fact the case. Remarkably, the same effect is directly observed in the gluon Wigner and Husimi distributions computed in Sec. III. We further examine and validate this interpretation by studying variations of the proton size in Appendix C.
Finally, we study the dependence of the dijet cross section and v 2 on the photon virtuality Q 2 . In Fig. 18 we show longitudinal and transverse cross sections as a function of Q 2 , where jPj ¼ 1 GeV, Δ ¼ 0.1 GeV and W ¼ 100 GeV, integrated over z and the azimuthal angles, θðΔÞ and θðPÞ. We choose a CGC parametrization with m ¼ 0.2 GeV,m ¼ 0.4 GeV, c ¼ 0.75 and α s ¼ 0.21 (cf. black curves in Figs. [15][16]. While the longitudinal cross section vanishes at Q 2 ∼ 0 (photoproduction), it becomes dominant for Q 2 ≥ 30 GeV 2 . Overall, the total cross section decreases at large Q 2 . The inset shows the ratio of the longitudinal to the transverse cross section.   In Fig. 19, we show the Q 2 dependence of the transverse, longitudinal and total v 2 with kinematics as in Fig. 18. Individually, both the transverse and longitudinal v 2 decrease with Q 2 . This is understood, since small dipoles contribute at large Q 2 which are less sensitive to correlations in the gluon dipole distribution. Interestingly, the v 2 of the transverse component turns negative at approximately Q 2 ¼ 15 GeV 2 .

V. CONCLUSIONS
We have computed gluon Wigner and Husimi distributions of the proton from the color glass condensate effective theory, at leading logarithmic order and including energy evolution by means of the JIMWLK renormalization group equations. We have studied angular correlations between the impact parameter and transverse momentum in Wigner and Husimi distributions, as well as between the impact parameter and dipole orientation in the dipole amplitude at small x. Our results qualitatively agree with the model calculations of Ref. [23], but differ at small jPj where nonperturbative contributions are most important. We studied the energy dependence of azimuthal correlations in gluon distributions, finding a decrease of these correlations with decreasing x, related to the geometric growth of the proton. We investigated differences between Wigner and Husimi distributions and pointed out difficulties for a physical interpretation of the latter due to its dependence on a smearing parameter.
A possible experimental consequence of elliptic correlations in the gluon Wigner distribution is corresponding elliptic modulations in the diffractive dijet production cross sections in e þ p collisions. We computed these within the color glass condensate framework for typical electron-ion collider energies in the correlation limit and studied angular correlations between the dijet transverse momentum and target recoil.
Starting from a simple modification of the IP-Sat model, we first established a baseline to disentangle purely kinematical correlations from those originating in the gluon distribution. We found a sizable hcos θðΔ; PÞi modulation in addition to the expected elliptic contribution. This kinematic effect is due to the energy dependence of the dipole amplitude and vanishes asymptotically in the correlation limit. Further, we presented dijet cross sections for charm jets, avoiding nonperturbative contributions from large dipoles.
We computed charm-dijet cross sections directly from the CGC by solving the JIMWLK renormalization equations for fundamental Wilson lines with initial conditions constrained by HERA data at relatively large x P . We studied the dependence of the elliptic correlation on relative transverse (jPj) dijet momentum, separately for both photon polarizations. We predicted the energy dependence of the elliptic anisotropy parameter and found a decrease towards smaller x, related to the growth of the proton with energy.
We computed the Q 2 dependence of the dijet cross section and elliptic modulation and observed a decrease of longitudinal and transverse elliptic modulations with photon virtuality, due to the suppression of larger dipoles in the photon wave function with increasing Q 2 . The transverse v 2;T turns negative at approximately Establishing a direct quantitative equivalence between correlations in dijet cross sections and gluon Wigner distributions is possible in restricted cases, such as for dijet production from real (Q 2 ∼ 0) photons in p þ A collisions [97]. In contrast, at the electron-ion collider one probes a wide range of photon virtualities, beam energies and dijet kinematics. Our study illustrates that many kinematic effects must be considered when trying to extract information about gluon Wigner and Husimi distributions in the proton from dijet cross sections. This might be further complicated in a more refined phenomenological analysis, where fragmentation, detector acceptances and other effects are included, and by parameter dependences of those distributions themselves, as we illustrated for the Husimi distribution. Nevertheless, in our work we have predicted features of the elliptic modulation of gluon distribution functions, such as its energy dependence, which are manifest in the diffractive dijet cross section and can be easily verified in experiments. The extension of the presented formalism to nuclear targets is straightforward.
Constraining the gluon Wigner distribution at small x is also important for future experimental investigation of the proton spin structure at the EIC [10,[20][21][22]41]. Work to extend the CGC framework to consistently include parton spin in a generalized Wigner function formulation is ongoing [150,151].

ACKNOWLEDGMENTS
where z ¼ ðz − ; zÞ. Here, are Wilson lines in the fundamental representation.
Assuming xP þ ≈ 0 in the small-x regime and using the cyclicity of the trace, we write We now define z ¼ x − y and introduce an irrelevant center Further, using and replacing the quantum-mechanical matrix elements by a stochastic average in the CGC effective theory, where α ¼ g 2 =4π. A similar derivation was presented in Ref. [103]. An effective method to compute the Wigner distribution numerically was developed in Ref. [23]. We include only the first two nonzero harmonic components in the expansion, and write the Wigner distribution as In our numerical study we use α s ¼ 0.3.

APPENDIX B: COORDINATE CONVENTIONS
In this Appendix, we discuss different transverse coordinate definitions to study angular correlations in coherent diffractive dijet production. In the main part of this manuscript we follow the convention of Ref. [7] [Eq. (9)], but a different convention was used in Ref. [132] (albeit for a different process),P ≡zp 0 − zp 1 ; For our study, in view of the kinematic effects discussed in Sec. IVA, an advantage of this convention is that the energy dependence of the dipole amplitude is independent of θðΔ;PÞ in these coordinates, One may conclude that these variables are more suitable for our study, potentially allowing to extract the gluon distribution at a fixed x P . However, now the cross section for dijet production cannot be written as a Fourier integral over the product of the photon wave function and dipole distribution. This is because r, b andP,Δ are not Fourier conjugates. Instead, if we naively evaluate the cross sections (10) and (12) for the IP-Sat model at fixed jPj, jΔj we find a nonzero v 2 even if there are no such correlations on the level of the gluon dipole distribution. This result is illustrated in Fig. 20, where we plot the normalized dijet cross section at jΔj ¼ 0.1 GeV, jPj ¼ 1 GeV, integrated over z ∈ ½0.1; 0.9 and with photon and proton beam kinematics as in Fig. 8. Figure 21 shows the z dependence of the v 2 using this convention, where nonzero dijet correlations are seen even in the absence of correlations in the gluon dipole distribution. We assume no angular correlations between the impact parameter and dipole orientation, yet a nonzero v 2 is seen.
To remedy this situation, one might alternatively parametrize the gluon Wigner distribution in Fourier conjugates to Eq. (B1),b thereby mixing longitudinal and transverse coordinates.

APPENDIX C: PROTON SIZE DEPENDENCE
An important consequence of JIMWLK evolution is the growth of the proton with energy. In this Appendix we investigate the dependence of v 2 on proton size. Here, instead of evolving using the JIMWLK equations, we directly match the x P dependence of the saturation scale Q s from the IP-Sat model, cf. Eqs. (19)- (20). In this matching procedure, we manually change the proton profile (16) by hand thus mimicking a geometric growth of the target, while keeping the saturation scale fixed.
Our results for jPj ¼ 1 GeV, jΔj ¼ 0.1 GeV, W ¼ 100 GeV, Q 2 ¼ 1 GeV 2 are shown in Fig. 22. Here we plot the z-integrated total v 2 as a function of the proton size in Eq. (16) for B p ¼ 4, 6, 16 GeV −2 . For larger B p the gradients of the color charge density become smoother, thus reducing the dependence of the cross section on the relative orientation between the dipole orientation and impact parameter.
These results suggest that the growth of the proton with energy is the most important effect behind the dependence shown in Fig. 17, while the change of the saturation scale Q s with energy has a minor effect on v 2 .