The importance of kinematic twists and genuine saturation effects in dijet production at the Electron-Ion Collider

We compute the differential yield for quark anti-quark dijet production in high-energy electron-proton and electron-nucleus collisions at small $x$ as a function of the relative momentum $\boldsymbol{P}_\perp$ and momentum imbalance $\boldsymbol{k}_\perp$ of the dijet system for different photon virtualities $Q^2$, and study the elliptic and quadrangular anisotropies in the relative angle between $\boldsymbol{P}_\perp$ and $\boldsymbol{k}_\perp$. We review and extend the analysis in [1], which compared the results of the Color Glass Condensate (CGC) with those obtained using the transverse momentum dependent (TMD) framework. In particular, we include in our comparison the improved TMD (ITMD) framework, which resums kinematic power corrections of the ratio $k_\perp$ over the hard scale $Q_\perp$. By comparing ITMD and CGC results we are able to isolate genuine higher saturation contributions in the ratio $Q_s/Q_\perp$ which are resummed only in the CGC. These saturation contributions are in addition to those in the Weizs\"ackerWilliams gluon TMD that appear in powers of $Q_s/k_\perp$. We provide numerical estimates of these contributions for inclusive dijet production at the future Electron-Ion Collider, and identify kinematic windows where they can become relevant in the measurement of dijet and dihadron azimuthal correlations. We argue that such measurements will allow the detailed experimental study of both kinematic power corrections and genuine gluon saturation effects.


Introduction
In the past decades, high energy collider experiments have successfully verified that quantum chromodynamics (QCD) is the theory of strong interactions of quarks and gluons (partons) inside hadrons and nuclei. A systematic program to uncover the structure of hadrons began with the introduction of parton distributions functions (PDFs) which describe the parton densities as functions of the longitudinal momentum fraction x. Beyond this onedimensional picture, transverse momentum dependent (TMD) PDFs have been introduced to characterize the three-momenta of partons inside hadrons [2][3][4]. Complementarily, generalized parton distributions (GPDs) have been defined by furnishing PDFs with the two dimensional transverse spatial distribution of partons resulting in a tomographic picture of hadrons and nuclei [5][6][7].
Among these various processes, forward particle azimuthal angle correlations are powerful observables to access the small-x structure of hadrons and nuclei at current and future collider experiments [60]. A paradigmatic example is that of inclusive dihadron production in d + Au collisions at RHIC, where a suppression in the back-to-back peak relative to p + p collisions [61,62] might signal the presence of gluon saturation [35,63,64] 1 . At the future Electron-Ion Collider (EIC) [66][67][68] and at the LHeC/FCC-he [69] a similar measurement has been proposed where a depletion of the back-to-back peak in forward (electron-going) dihadron production is expected to be observed when going from proton to nuclear deeply inelastic scattering (DIS) [37,70]. This measurement at the EIC offers two advantages: i) there is no need to assume a hybrid dilute-dense factorization in which one convolutes with PDFs of the proton since in DIS the projectile is a virtual photon whose kinematics can be reconstructed by measuring the scattered electron, ii) in the correlation (back-toback kinematics) limit, it is only sensitive to one type of gluon TMD distribution 2 : the Weizsäcker-Williams (WW) type. Furthermore, measurements of the azimuthal distribu- 1 Another mechanism that could provide a suppression of the back-to-back peak is the momentum broadening due to cold nuclear matter energy loss and coherent power corrections, which has been studied for dihadron production in p + p and d + Au in [65]. 2 Unlike p/d + p/Au, which is sensitive to eight different types of gluon TMDs [64,70,71].
tion of the momentum imbalance k ⊥ of the dihadron/dijet system provide the opportunity to access the linearly polarized WW gluon TMD. In particular, the elliptic anisotropy in the angle between the momentum imbalance k ⊥ and the relative dijet momentum P ⊥ is proportional to the ratio of linearly polarized to unpolarized gluon pairs [72][73][74]. Both unpolarized and linearly polarized WW gluon TMD are expected to be sensitive to saturation effects at small x as they resum contributions in powers of Q s /k ⊥ , making the studies of dihadron and dijet anisotropies promising channels to investigate gluon saturation [37,38]. While the TMD framework for inclusive dijet production is expected to hold near backto-back kinematics, there are important higher twist corrections that must be resummed for more controlled phenomenological predictions. At small x, the CGC EFT approach to highenergy QCD [18][19][20]22] encodes these higher order corrections in the dipole and quadrupole correlators of light-like Wilson lines 3 . In [1], we computed the inclusive dijet production differential cross section in the CGC formalism employing the Gaussian approximation [70,75,76] and at realistic energies for the future EIC, and observed deviations from the TMD framework at large momentum imbalance k ⊥ ∼ P ⊥ . In addition, we also found differences at small momentum imbalance if also Q s ∼ P ⊥ , which were enhanced in nuclear DIS due to the enhanced nuclear saturation scale Q 2 s ∝ A 1/3 . These genuine saturation contributions are different than those existing in the WW gluon TMD at small x which scale as Q s /k ⊥ [70,74].
When studying power corrections to semi-inclusive processes with several hard and semi-hard scales involved, it is necessary to distinguish terms according to the ratios they resum. Indeed, the presence of a semi-hard k ⊥ and a saturation scale Q s in addition to the hard scale Q ⊥ implies that power-suppressed terms can scale as k ⊥ /Q ⊥ or Q s /Q ⊥ , which can be significantly enhanced when compared to the simpler Λ QCD /Q ⊥ twist corrections. This approach has been recently applied to the production of quark anti-quark dijets in p+p and p + A collisions in [77]. The present manuscript follows the same spirit of [77] combined with the results in [78,79], which allow us to make the distinction of higher twists into kinematic (k ⊥ /Q ⊥ ) and genuine (Q s /Q ⊥ ) saturation contributions for the electroproduction of a quark anti-quark dijet. On the level of the differential cross section they can be schematically identified as follows: where the hard scale Q ⊥ of the process is built from the relative momentum P ⊥ and the virtuality of the photon Q, both controlling the dipole size r ⊥ . For simplicity, one can define Q ⊥ = max(Q, P ⊥ ).
In practice, the so-called genuine saturation terms correspond to the contribution of genuine higher twist operators. At small x, they correspond to operators with higher powers of the gluon field strength tensor in the form of gF µν insertions along with gauge links that maintain gauge invariance. As we will prove in this article, these corrections yield powers of Q s /Q ⊥ in the CGC because of semi-perturbative gluon loops whose internal transverse momenta peak around Q s . The ITMD framework resums only the kinematic power corrections (kinematic twists) leaving aside the genuine saturation contributions. This framework was originally introduced to study dijet production in proton-nucleus collisions [71,80,81], and recently extended for DIS and photo-production ( [78,79]). By comparing TMD, ITMD and CGC frameworks, we are able to identify the separate roles of kinematic and genuine higher twists in the production of quark anti-quark dijets at the EIC, which is the goal of the present manuscript.
This manuscript is organized as follows. In Sec. 2 we begin by briefly reviewing the results in [82] relating the product of two light-like Wilson lines to a transverse gauge link. This allow us to generalize the small dipole size expansion in [70], and isolate kinematic and genuine higher twists. In Sec. 3, we discuss the inclusive production of dijets in high energy (small-x) DIS. We start with the computation in the CGC formalism in momentum space, where we also take the opportunity to set up the notations. We choose to express the differential cross-section in terms of the momentum imbalance k ⊥ and the relative momentum P ⊥ of the dijet pair. Then, we review the derivation of the TMD limit for back-to-back dijets (k ⊥ , Q s Q ⊥ ) and introduce the WW gluon TMD. Next, we study the ITMD framework which resums the kinematic twists (in powers of k ⊥ /Q ⊥ ), and obtain analytic expressions for the hard factors. We briefly discuss genuine higher twists and argue that these are parametrically of the order Q s /Q ⊥ . In Sec. 4, we discuss the setup for our numerical computation: the initial conditions, the small-x evolution, the Gaussian approximation for high energy CGC correlators and the WW gluon TMD, and the computation of the harmonics of the differential cross-section in the angle φ = φ P ⊥ − φ k ⊥ . Our numerical results for the production of quark anti-quark dijets at realistic EIC kinematics are presented in Sec. 5. We focus on the hadronic sub-amplitudes γ * L/T + p/Au → qq + X, and compare the results in the TMD, ITMD and CGC formalisms. We first present the angle averaged yield (over φ P ⊥ and φ k ⊥ ) as a function of momentum imbalance k ⊥ at different values of P ⊥ and virtuality Q. Then, we focus on the production close to back-to-back kinematics at various values of virtuality and relative jet momentum. Finally, we present our results for the elliptic and quadrangular azimuthal anisotropies in the angle φ. In Sec. 6 we summarize our main results and comment about their implications to the measurements of azimuthal dijets and dihadron correlations at the EIC.

Kinematic and genuine higher twists in the Wilson line pairs
In the CGC, the scattering of a colored particle moving along the minus light-cone direction through the background field A µ of a nucleus propagating in the plus light-cone direction can be characterized (in the eikonal approximation) as a color rotation of the particle given by the light-like Wilson line: where t a are the generators of SU (3) in the fundamental representation, and A +,a is in Lorenz gauge ∂ µ A µ = 0. Here P stands for path ordering such that the operator at z = −∞ is in the rightmost position, while that at z = +∞ is in the leftmost position 4 . For the sake of clarity, we will decorate the gauge field with a tildeÃ when working inÃ + = 0 light-cone gauge as opposed to Lorenz gauge ∂ µ A µ = 0. At the level of the amplitude, the inclusive production of a quark anti-quark dijet in the CGC will contain the product of two Wilson lines V (x ⊥ )V † (y ⊥ ) corresponding to the color rotation of the quark and the anti-quark, respectively. The seminal work in [70] established the connection between CGC and TMD amplitudes (cross-sections) by noting that the first term in the expansion in powers of r ⊥ = x ⊥ − y ⊥ of this product produces the transverse gauge fieldÃ i ⊥ in light-cone gaugeÃ + = 0: 3.1 and 3.2). Such transverse gauge fields can be written in terms of the field strength tensor F i+ , which constitute the building blocks of small-x gluon TMDs (see Appendix F).
It seems natural to wonder what is the behavior of the higher terms in powers of r ⊥ in the expansion in Eq. (2.2). This is the approach followed in [83], where part of the quadratic term in r ⊥ has been isolated. However, one quickly realizes that order by order in the r ⊥ expansion, the resulting terms are cumbersome to organize in a way which preserves explicit QCD gauge invariance for the involved operators. Recent developments on the CGC/TMD correspondence rely on the possibility to express the product of two Wilson lines as a transverse gauge link at x − = −∞ [82]: One can then show (see Appendix A for details) that as a consequence of recursive relations satisfied by the transverse gauge link, it is possible to re-organize all terms beyond the linear expansion in gÃ i ⊥ as follows: . (2.4) We note that this is an exact relation, and that it is organized in a gauge invariant way circumventing the difficulty introduced by the naive expansion in the dipole size r ⊥ . The infinite series obtained via iteration of the right-hand side of Eq. (2.4) re-expresses the multiple scattering of a quark anti-quark pair as a combination of gluon saturation (embodied in the strong transverse gauge fieldÃ i ⊥ (z ⊥ ) ∼ 1/g) and multiple scattering in light-cone gauge 5 . In Sec. 3.3 we will show explicitly that when the linear term in Eq. (2.4) is used to calculate the amplitude for quark anti-quark dijet production, it will only involve the WW gluon TMD (single scattering) but with a modified hard factor that resums all kinematic power corrections in k ⊥ /Q ⊥ (ITMD). On the other hand, the higher order terms in the quadratic term (gÃ ⊥ ) 2 (higher order operators in gF i+ ) resum contributions in powers of Q s /Q ⊥ , corresponding to double scattering (and beyond through iteration), which will be discussed in Sec. 3.4.

Inclusive dijet production in high energy DIS
In this section, we review the computation of the differential cross section for the production of a forward quark anti-quark dijet in DIS within the CGC EFT e(k e ) + A(P A ) → e(k e ) + q(k 1 ) +q(k 2 ) + X, (3.1) where A can be proton or a nucleus. The kinematic variables are detailed in Table 1.
The leptonic and hadronic parts of this process can be separated by expressing the production cross-section as follows: where λ is the polarization of the exchanged virtual photon, Q 2 is its virtuality, and W 2 is the center of mass energy per nucleon of the γ * − A system. The leptonic part is contained in the photon fluxes: where the inelasticity is given by y = W 2 +Q 2 −m 2 n s−m 2 n and √ s is the center of mass energy per nucleon in the e + A collision. 5 We note that the physical picture of gluon saturation is gauge dependent. In Lorenz gauge ∂µA µ = 0, high gluon-density manifests as multiple scattering encoded in the Wilson lines. On the other hand, in lightcone gaugeÃ + = 0 the phenomenon of gluon saturation becomes manifest in the WW gluon distribution, which represents the gluon number density (see e.g. Sects. 2.5 and 2.6 in the review article [18]). incoming (outgoing) electron four-momentum q = k e − k e virtual photon four-momentum k 1,2 quark (anti-quark) four-momentum z 1,2 quark (anti-quark) longitudinal momentum fraction η 1,2 quark (anti-quark) rapidity quark anti-quark dijet transverse momentum imbalance P ⊥ = z 2 k 1⊥ − z 1 k 2⊥ quark anti-quark dijet relative transverse momentum s = (P n + k e ) 2 nucleon-electron system center of momentum energy squared W 2 = (P n + q) 2 nucleon-virtual photon system center of momentum energy squared m 2 n = P 2 n nucleon invariant mass squared M 2 qq = (k 1 + k 2 ) 2 invariant mass squared of the dijet system.
virtuality squared of the exchanged photon In the present work, we focus on the hadronic content of the scattering encoded in the sub-process: We will work in light-cone coordinates 6 , and in a frame where the proton or nucleus propagates in the plus light-cone direction and the virtual photon in the minus light-cone direction: with m A the proton (nucleus) mass. The momenta of the produced q andq are given by respectively. We will denote the longitudinal momentum fractions z 1,2 = k − 1,2 /q − , and we will neglect the quark masses. 6 We follow the convention a ± = 1 √ 2 (a 0 ± a 3 ), such that a µ bµ = a + b − + a − b + − a ⊥ · b ⊥ , where a ⊥ and b ⊥ are the 2-dimensional transverse vectors. We will denote the magnitude of transverse vectors as a ⊥ = (a ⊥ · a ⊥ ) 1/2 .

Inclusive quark anti-quark production in the CGC
In the CGC EFT, the large-x degrees of freedom in the nucleus are treated as stochastic color sources ρ a A which generate the strong classical gluon field A µ , characterizing the smallx gluon content of the nucleus.
For a fast moving nucleus along the plus light-cone direction, the color sources generate a current of the form which is independent of light-cone time x + . It can be easily verified that the Yang-Mills equations [D µ , F µν ] = J ν are satisfied by the gauge field where α(x − , x ⊥ ) solves the Poisson equation The solution in Eq. (3.11) is in Lorenz gauge ∂ µ A µ = 0. An infinite set of solutionsÃ can be found by application of gauge transformations One notable choice is the light-cone gaugeÃ + = 0 , which can be obtained from the gauge transformation In such a case one obtains a solution in terms of a purely transverse gauge field: Note that this choice relies implicitly on the boundary conditionÃ i ⊥ (+∞, x ⊥ ) = 0 which is allowed by the freedom of the residual subgauge condition in light-cone gauges. It will be useful to denote the transverse gauge field at x − = −∞ by simply suppressing the light-cone argument: The expectation value of any observable in the CGC is then computed in perturbation theory in the presence of the background field A µ for a given configuration of sources ρ A , and then averaging over all possible configurations: where W Y [ρ A ] is a gauge invariant weight functional for the configuration of ρ A , at the rapidity scale Y = log(1/x) (see also Ref. [84] for a discussion of different evolution rapidity variables). Here x is the typical longitudinal momentum fraction probed by the observable. Following this prescription, the differential cross-section reads: where M λ CGC [ρ A ] is the CGC amplitude for the production of quark anti-quark in the collision of a virtual photon γ * with the external classical background field (shock-wave) produced by the source configuration ρ A characterizing the nucleus A.

Amplitude
At leading order in the CGC, the forward production of a dijet pair proceeds by the splitting of the virtual photon γ * into a quark anti-quark dipole that subsequently scatters off the field A µ , and then fragments into two jets (see Fig. 1  The result for the computation of quark anti-quark production in DIS within the CGC is well-known and has been obtained in [70,72]. Here we briefly review the computation of the amplitude for this process following momentum space Feynman rules [86]. The expression for the amplitude is where we followed the standard momentum space Feynman rules of QCD+QED in ∂ µ A µ = 0 gauge for QCD, supplemented by the effective CGC vertex for quark interaction with the back-ground field. The Feynman propagator for a free massless quark reads S 0 (l) = i / l l 2 + i , (3.20) 7 We expect that the contribution from the quark initiated channel γ * + q → q + g to be negligible in low-x kinematics because of the dominance of gluon exchanges. For a recent study of this process and the impact of gluon saturation and multiple parton scattering see [85]. and the effective CGC vertex for the eikonal multiple scattering of a quark (anti-quark) off the color field A µ [87,88]: where i, j represent the color indices, and l and l are the outgoing and incoming momenta of the quark, and the light-like Wilson line in the fundamental representation was introduced in Eq. (2.1). The superscript sgn(p) in Eq. (3.21) denotes matrix exponentiation: where the latter follows from the unitarity of V (z ⊥ ).
It is convenient to work in A − QED = 0 light-cone gauge for the QED part of the amplitude and in a frame where the transverse momentum of the photon is 0 ⊥ . The polarization vectors for the virtual photon read 8 where ±1 ⊥ = 1 √ 2 (1, ±i), λ = 0 denotes the longitudinal polarization, and λ = ±1 the two transverse polarizations.
Subtracting the non-interacting piece and factoring out an overall delta function 2πδ(k − 1 + k − 2 − q − ), the amplitude can be organized as follows with the perturbative factor: where σ and σ are the helicities of the quark and anti-quark respectively. The computation of the perturbative factors is straightforward, we briefly outline it in Appendices C and D. We find where ε 2 = z 1 z 2 Q 2 . 8 In the Lorenz gauge ∂µA µ = 0 for QED, the longitudinally photon polarization takes the form λ=0,µ =

Differential Cross-section
The color structure obtained by squaring the amplitude, summing over the colors of the quark and the anti-quark in the final state, and averaging over the different color charge configurations ρ A is given by where we define the CGC average of the dipole and quadrupole operators respectively as The impact factor is obtained by squaring the perturbative factor in Eq. (3.25) and summing over the helicities of the quark and the anti-quark in the final state: The differential cross-section for quark anti-quark production in the CGC is then is an overall minus light-cone momentum conserving delta function, and for convenience we defined the measure: The impact factors corresponding to longitudinally and transversely polarized photons can be easily obtained by inserting Eqs. (3.27) and (3.28) in Eq. (3.32): For transversely polarized photon contributions we averaged over both polarizations λ = ±1. The obtained differential cross section agrees with the previous calculations presented in Ref. [70,72,73]. A convenient choice of momentum variables to characterize the production of dijets is 38) which are conjugate to the coordinate variables: respectively the dipole impact parameter and the dipole relative vector. The momentum imbalance k ⊥ measures the deviations of the dijet system from being back-to-back in the transverse plane, while the relative momentum P ⊥ is closely related to the invariant mass of the dijet pair This pair of momentum variables is also useful to establish the relation between the CGC EFT and the TMD framework as will be shown in the next section.
With this alternative choice of coordinates, the differential cross-section is written as where the measure is defined as with the right-hand side defined in Eq. (3.29).

The TMD limit
The transverse momentum dependent (TMD) framework for the inclusive production of a quark anti-quark pair is expected to hold when the momentum imbalance is small relative to the typical transverse momenta of the q andq, i.e., k ⊥ P ⊥ . In coordinate space, this condition is expected to be equivalent to the small dipole expansion r ⊥ b ⊥ (see [70]), although this assertion will be corrected in Section 3.4. Such an expansion applied to the Wilson line correlators in the CGC amplitude (Eq. (3.24)) results in: The correlator appearing on the right hand side is the pure gauge transverse gauge field one would obtain if working in theÃ + = 0 gauge [89]: as introduced in Sec. 3.1.
Inserting the result of Eq. (3.45) in Eq. (3.24) and ignoring higher order terms in r ⊥ (or equivalently in 1/Q ⊥ ), we obtain the amplitude in the TMD framework: The perturbative factors in Eq.(3.48) for longitudinally and transversely polarized photons are given respectively by: Then the differential cross-section reads: where the hard factor is 53) and the soft factor is the Weizsäcker-Williams [70,90] gluon TMD defined as (see Appendix F for more details): (3.54) While at large k ⊥ , one expects a power tail xG ii (x, 72]. This behavior at low momentum imbalance results in a nuclear suppression of back-to-back dihadron production in DIS [37].
The differential cross-section for quark anti-quark production in DIS within the TMD limit is given by [70,72,73]: where φ is the angle between P ⊥ and k ⊥ , and the WW correlator has been decomposed into trace and traceless parts: The trace and traceless parts are known as unpolarized and linearly polarized gluon WW TMD distributions, respectively. They obey the inequality which renders the cross-sections in Eqs. (3.55) and (3.56) positive [3]. Note that in the TMD framework the elliptic anisotropy defined as is proportional to the ratio of the linearly polarized to the unpolarized distribution [72,74]:

The improved TMD framework: resumming kinematic twists
While the CGC EFT takes into account all possible twist corrections that are not suppressed by a factor of 1/W , the so-called small-x improved TMD framework (ITMD for short) [80] only resums specific power corrections which correspond to kinematic corrections in the hard sub-amplitude: while the TMD framework expresses the observable as the ITMD framework rewrites it as . It is crucial that the distribution xG ij is the same for both schemes. Indeed, the ITMD framework neglects higher order corrections gÃ i ⊥ in the operators that appear in the CGC EFT. Unlike the TMD approximation, where one expands in r ⊥ , our starting point is the truncation of the expansion introduced in Sec. 2 (see Eq.(2.4)): This expression can be cast into the following form (see Appendix A): (3.65) When plugged into the amplitude from Eq. (3.24), it yields: whereÃ i ⊥ (k ⊥ ) has been defined in Eq. (3.49) and the hard factor is defined as Note that by neglecting powers of k ⊥ · r ⊥ (or equivalently powers in k ⊥ /Q ⊥ ) one recovers Eq. (3.48) and thus Eq. (3.66) reduces to the TMD limit in Eq. (3.47). The differential cross-section for the production of an quark anti-quark pair in the ITMD framework reads: The computation is almost identical to that in the TMD framework, except for the hard factor in Eq. (3.69), whose dependence on k ⊥ resums kinematic power corrections in k ⊥ /Q ⊥ . The hard factor in the ITMD framework (Eq. (3.67)) can be computed fully analytically (see Appendix E for details). When the virtual photon is longitudinally polarized, we find , (3.70) and for the transversely polarized photon, we have where we used k 1⊥ = P ⊥ + z 1 k ⊥ , k 2⊥ = −P ⊥ + z 2 k ⊥ (see Eqs. (3.37) and (3.38)), and we defined: These expressions will be the ones used in our numerical evaluation of the ITMD differential cross-section 9 . Some interesting limits are discussed in Appendix E.

Genuine higher twists
In the landmark article on the equivalence between CGC observables and TMD-factorized observables [70], the notion of the correlation limit was introduced as a means to justify taking the Taylor expansion from which TMD distributions emerge as correlators of Wilson lines and derivatives of Wilson lines. While this limit does correspond to the TMD limit, there is an additional subtlety worth discussing. Indeed, it is common to identify the correlation limit with the back-to-back kinematic limit. Given that r ⊥ and b ⊥ are respectively Fourier conjugated to P ⊥ and k ⊥ , in the k ⊥ P ⊥ limit one can justify assuming r ⊥ b ⊥ . More recent insights on the CGC/TMD correspondence [79] have revealed the limit of this hypothesis. Because in CGC observables dipole sizes appear both in the non-perturbative correlators and in the hard subamplitudes, powers of r ⊥ can actually be enhanced by powers of the saturation scale Q s via non-perturbative effects in the target. Assuming that the k ⊥ P ⊥ limit (or more precisely k ⊥ Q ⊥ ) and the r ⊥ b ⊥ limit are indistinguishable is tantamount to neglecting all powers of Q s /Q ⊥ . As we will show in this section, there is actually a non-zero contribution from genuine higher twists even in the k ⊥ → 0 limit. These Q s /Q ⊥ corrections are taken into account neither in the TMD nor ITMD framework. 10 The quadratic term in Eq. (2.4) can be cast into the following form (see Appendix A): The last line in the r.h.s. of this equation is a building block to construct genuine higher twist TMD distributions, e.g. gauge invariant distributions with 4 field strength tensors carrying non-zero transverse momenta, but it is worth noting that the V (v 1⊥ )V † (v 2⊥ ) operator could also be recursively expanded in powers of gÃ i ⊥ ad nauseam. When plugged into the 9 While this work was in progress, the hard ITMD factors were computed independently in [91]. We have verified that both results are equivalent. 10 Recall the saturation scale Qs appears in the TMD (and ITMD) framework only as powers Qs/k ⊥ in the WW gluon TMD.
CGC amplitude from Eq. (3.24), we obtain the following (g 2Ã2 ⊥ )-suppressed contribution: We can now see the mechanism through which a Q s /P ⊥ correction emerges from such a contribution. The 2 gluons in the amplitude both carry transverse momenta, whose sum is k ⊥ . The momentum difference ⊥ , which is Fourier conjugated to the difference in positions v 1⊥ and v 2⊥ , provides an additional loop variable. In CGC descriptions for the target, the momentum integral will typically peak around | ⊥ | ∼ Q s , while due to the hard factors the dipole sizes peak around r ⊥ ∼ 1/Q ⊥ . Overall, every power of ( ⊥ · r ⊥ ) yields a Q s /Q ⊥ correction, regardless of the magnitude of k ⊥ . Taking the limit k ⊥ → 0 ⊥ (back-to-back) in the hard subfactor, we get the non-zero correction We finally see the difference between the so-called correlation (small dipole) limit and the kinematic back-to-back limit: genuine higher twists contribute in the latter but not in the former. Let us illustrate this more concretely by taking the leading power of Q ⊥ in Eq. (3.75). From the previous discussion, we know that each power of ( ⊥ · r ⊥ ) yields a power correction. From the leading term of such an expansion, we get It is particularly revealing to compare this to the TMD amplitude, given by We immediately see that in M λ g.h.t. , the hard factor is enhanced by a power of r ⊥ hence suppressed by a power of the hard scale after the integral is taken. This additional power is compensated by the fact that the operator now has its dimension (hence its twist) increased by 1 due to the presence of an additional gA j insertion. The action of a higher twist operator on target states will be one higher power of the typical target scale: Q s in the CGC or Λ QCD in dilute schemes.
To sum up this section: genuine higher twist corrections are suppressed by a power of the ratio between a target scale and a hard scale, and they contribute to the observable regardless of how k ⊥ compares to P ⊥ or Q. In the CGC, we expect the ratio to be linear in Q s , which means comparing CGC predictions to TMD predictions for dense targets in the exact back-to-back limit would provide a probe for genuine saturation effects.

Initial conditions and small-x evolution
The CGC target average of the dipole and quadrupole operators, Eqs. (3.30) and (3.31), are necessary ingredients to evaluate the dijet production cross section (3.33), and are defined as correlators of Wilson lines V (x ⊥ ). The Weizsäcker-Williams gluon distribution, Eq. (3.54), can also be written in terms of the dipole operator (and its derivatives) as we will show explicitly below. The energy (or momentum fraction x) dependence of the Wilson lines is given by the JIMWLK evolution equation [15-17, 89, 92-94]. In the so-called planar limit where the number of colors N c is considered large enough to neglect 1/(N 2 c − 1) corrections, it reduces to the Balitsky-Kovchegov (BK) evolution equation [95,96], which is a closed equation to describe the energy evolution of the dipole operator.
The JIMWLK evolution of Wilson lines on a transverse lattice can be solved numerically [97][98][99][100], which would allow one to directly evaluate the 2-and 4-point correlators, Eqs. (3.30) and (3.31). However, this is computationally demanding, and instead we will employ the Gaussian approximation discussed in Sec. 4.2 which allows one to express higher point correlators in terms of the dipole correlator only 11 .
We use the dipole-target scattering amplitude obtained in Ref. [32] (see also Ref. [23]). Here, at the initial x = 0.01 the functional form of the dipole operator (in case of proton targets) is obtained from a McLerran-Venugopalan model [12] based parametrization and written as (see also Refs. [101,102]) The dipole amplitude at x < 0.01 is obtained by solving the BK equation over Y = ln 0.01 x units of rapidity, including running coupling corrections [103]. The free parameters (Q 2 s,0 , the scale of the running coupling in coordinate space and the proton transverse area S ⊥ ) are determined in Ref. [32] by performing a fit to the proton structure function data from HERA [104]. We note that there has recently been progress to promote the structure function calculations to next-to-leading order accuracy [24,[105][106][107][108][109][110][111], but for consistency we use the leading order result from Ref. [32] in our leading order calculations.
To generalize the dipole-proton operator S (2) to describe the interaction with a heavy nucleus, we again follow Ref. [32]. At the initial condition, the dipole operator is written as with S ⊥ = 18.81 mb , and T A is the spatial density profile obtained from the Woods-Saxon distribution. To obtain the dipole amplitude at smaller x, the BK equation is solved at fixed impact parameter. This approach results in vanishing nuclear effects in the dilute region and successful phenomenology, see e.g. Refs. [42,112,113]. Since the determination of the impact parameter is subtle in electron-nucleus collisions (see however Refs. [114,115]), in this work (as in Ref. [1]) we use an impact parameter independent dipole amplitude, obtained by evaluating the dipole at the median impact parameter , as a proxy for minimum bias collisions. This results in a nuclear oohmp factor of the effective initial nuclear saturation scale relative to the proton saturation scale. For gold nuclei, we find S ⊥ AT A ( b ⊥ ) ≈ 3.1 .

Gaussian approximation for high energy correlators
Assuming that the color charges in the target are Gaussian distributed both at the initial condition and after the small-x evolution, it becomes possible to express all higher point correlators in terms of the two-point correlator only. This is referred to as a Gaussian approximation or Gaussian truncation [70,75,76], and is used in this work to express the quadrupole operator (3.31) in terms of the dipole operator (3.30) satisfying the BK evolution as discussed in the previous Section. The validity of the Gaussian approximation has been numerically confirmed in Ref. [116] by comparing the approximatively calculated quadrupole operator to the one obtained from JIMWLK-evolved Wilson lines, and analytically justified in [117]. It has also been used for example in phenomenological analyses of multi particle production [35,118,119] and to evaluate the higher point correlators in the next-to-leading order BK evolution equation [110,120] in Ref. [121]. Following Ref. [70], the quadrupole operator in the Gaussian approximation can be written as Here we used the following definitions: where C F = (N 2 c − 1)/(2N c ) = 4/3 is the fundamental Casimir. Note that we assumed that the dipole only depends on its relative vector and not the impact parameter.
The Weizsäcker-Williams distribution from Eq. (3.54) can be conveniently written in terms of Γ and its derivatives as [70,83,122] where C A = N c = 3 is the adjoint Casimir, and S ⊥ represents the transverse area of the target, which factors out when the dipole is translationally invariant (impact parameter independent). The trace and traceless components xG 0 and xh 0 defined in Eq. (3.57) can now be written as With these results, it becomes possible to evaluate all cross sections in terms of the BK evolved dipole operator S (2) only.

Computing harmonics
Before we proceed, we note that due to the translational invariance of the dipole, the differential cross-section is proportional to the overall area of the target (proton/nucleus) S ⊥ . Thus, we shall study the differential yield Due to overall rotational invariance, the differential yield is independent of the angle Φ = φ P ⊥ + φ k ⊥ ; thus it is sufficient to characterize the transverse momenta of the jets with P ⊥ , k ⊥ and the relative angle φ = φ P ⊥ −φ k ⊥ . It is then convenient to decompose the differential yield in angular modes with respect to φ : where the modes are given by From these quantities, we can then compute the elliptic and quadrangular anisotropies: In the TMD limit the only non-vanishing mode is v 2,λ for which there are explicit expressions in terms of the kinematic variables and the WW gluon TMD (see Eqs. (3.60) and (3.61)).
In the ITMD and CGC there are no simple expressions for the anisotropies, thus they have to be computed numerically by first evaluating the differential yield as a function of relative angle φ, and subsequent integration. While this is the approach we follow to compute these modes in the ITMD framework, for the computation in the CGC it is advantageous to use the following identity: Inserting the CGC differential cross-section Eq. (3.42) into Eq.(4.13) and using the identity above, we can perform analytically the integrals over φ P ⊥ and φ k ⊥ resulting in The computation of the modes and anisotropies for dijet production in the CGC reduces to the 6 dimensional integration above. Our results are shown in the next section.

Numerical results
In this section we numerically evaluate the differential yield for the inclusive quark antiquark production in proton and nuclear DIS. As argued in Sec. 4.3 (see Eq. (4.12)) we can write the differential yield as We will focus our study on the averaged angle differential yield N λ 0 (P ⊥ , k ⊥ ), the elliptic anisotropy v 2,λ (P ⊥ , k ⊥ ) and the quadrangular anisotropy v 4,λ (P ⊥ , k ⊥ ). The small-x BK evolution will be carried up to Y = log(0.01/x g ), where is the fraction of the target plus momentum transferred to the dijet system in the t-channel exchange. The left panel in Fig. 2 shows the saturation scale Q s as a function of x g for proton, and for gold at median impact parameter (as a proxy of minimum bias collisions) following our initial conditions (see Sec. 4.1) 12 . The saturation scale has been defined as where Y = log x 0 xg . The right panel in Fig. 2 displays the value of x g as a function of the kinematic variables. At the projected top EIC energies, requiring x g ≤ 10 −2 constrains the transverse momenta and virtualities of the dijet system as shown.
Our results will be shown separately for transversely and longitudinally polarized photons in collisions off protons and gold nuclei at a center of mass energy W = 90 GeV. We choose a configuration for the quark and the anti-quark in which they have identical longitudinal momenta z 1 = z 2 = 1 2 (note that the momentum fractions are related to pseudorapidities as z i = 2E n |k i⊥ |e −η i /W 2 where E n is the nucleon energy). Experimentally it is not directly possible to determine the photon polarization in dijet production, and as such the cross section can not be directly measured for the longitudinal and transverse 12 The modest values of the saturation scale displayed in Fig. 2 are a result of the parametrization of the dipole amplitude with MV initial conditions in [32]. Other parametrizations can result in larger values of the saturation scale. polarization states separately. As the photon flux factors f λ in the electron-nucleus cross section (Eq. (3.2)) depend on inelasticity, some insight can in principle be obtained by doing the same measurement at different √ s resulting in different inelasticities y. Additionally, as discussed in Ref. [38], the polarization dependent cross sections and elliptic modulations can be extracted from the data assuming that the functional form of the v 2,λ coefficients is know from theory. We will perform our computations in the TMD framework (Sec. 3.2), the improved TMD framework (Sec. 3.3), and the CGC EFT (Sec. 3.1). We remind the reader that by comparing the TMD and the improved TMD framework we gain access to kinematic power corrections: while the comparison of the CGC with the ITMD will help us assess the role of genuine saturation contributions: These saturation contributions are present in addition to those appearing in the WW gluon TMD that resum powers of Q s /k ⊥ .

Angle averaged differential yield
We start by presenting our results for the angle averaged differential yield. In Fig. 3 we show our results for N T 0 (P ⊥ , k ⊥ ) as a function of momentum imbalance k ⊥ , and at fixed virtuality Q 2 = 10 GeV 2 and relative momentum P ⊥ = 2 GeV.
In both proton and gold collisions, the CGC results in a suppression of ∼ 20% relative to the TMD when P ⊥ ∼ k ⊥ . This suppression is enhanced as the momentum imbalance k ⊥ is increased (away from the quark anti-quark back-to-back configuration) signaling the breakdown of the TMD framework. When kinematic power corrections in k ⊥ /P ⊥ are resummed into the hard factors of the improved TMD framework, we observe an excellent agreement between the CGC and ITMD results for k ⊥ > P ⊥ and up to large values of k ⊥ relative to P ⊥ , suggesting that genuine higher twist effects are suppressed for non-back-toback configurations and the suppression in the yield observed in the CGC relative to the TMD is driven by kinematic power corrections. Interestingly, the genuine saturation effects become visible in the back-to-back regime where deviations between TMD (or ITMD 13 ) and the CGC are observed, and which are enhanced from 7% in e + p to about ∼ 20% in e + Au collisions.
In Fig. 4, we present the ratios ITMD/CGC and TMD/CGC, and study their dependence as a function of k ⊥ at different virtualities Q and relative momenta P ⊥ . As expected, increasing either Q or P ⊥ eventually reduces the kinematic twists resulting in smaller differences between the CGC and TMD at a given k ⊥ . We also observe a systematic reduction of the genuine higher twist contributions at moderate P ⊥ and Q values. Thus we expect   that the improved TMD framework will provide good agreement with the CGC in the production of quark anti-quark pairs in DIS for kinematics in which P ⊥ or Q are significantly larger than the saturation scale. On the other hand, genuine twist corrections have a significant impact on the measurement of back-to-back dihadrons/dijets at low virtualities and transverse momenta in nuclear DIS where the saturation scale is enhanced. We now move to the differential yield N L 0 (P ⊥ , k ⊥ ) for the longitudinally polarized virtual photon case, shown in Fig. 5 as a function of momentum imbalance k ⊥ , and at a fixed virtuality Q 2 = 10 GeV 2 and relative momentum P ⊥ = 2 GeV. The results depicted are qualitatively similar to those in the transversely polarized case, but with some notable differences. We observe deviations between the CGC and the TMD at large momentum imbalance k ⊥ because of the importance of kinematic power corrections. Once these contributions are resummed in the ITMD the agreement with the CGC results is improved. We note that the agreement between CGC and ITMD at large momentum imbalance is not as good as in the transversely polarized case, which might be caused by the bias for larger dipole size contributions in the light-cone wave-function of longitudinally polarized photons 14 . This reveals that genuine saturation corrections also appear away from the back-to-back limit in the longitudinally polarized case. Ratios of ITMD and TMD to GCC for other choices of P ⊥ and Q 2 are shown in Fig. 6. We find both enhancement and suppression of the quark anti-quark production in the CGC relative to the ITMD, depending on the virtuality Q and the relative momenta P ⊥ . Overall we observe that the inclusion of kinematic power corrections improves the agreement with the CGC, and that increasing Q or P ⊥ reduces genuine higher twist contributions as one should expect. The apparent better agreement in Fig. 6 at P ⊥ = 3 GeV and Q 2 = 10 GeV 2 is coincidental as the CGC is transitioning from enhancement to suppression as Q 2 grows.
14 In contrast to fully inclusive DIS, in dijet production the longitudinal momentum fractions z1,2 are kept fixed. At the chosen values of z1 = z2 = 0.5 here (and at fixed virtuality Q 2 ), one can verify that the light-cone wave-function for a transverse photon ∼ K1(εr ⊥ ) grows more quickly than that of a longitudinal photon ∼ K0(εr ⊥ ) in the limit r ⊥ → 0. Figure 7. Ratio of TMD to CGC for near back-to-back kinematics k ⊥ ≈ 0 in γ * λ + Au → q +q + X, as a function of virtuality Q and relative momentum P ⊥ . Left: transversely polarized case. Right: longitudinally polarized case.
In order to single out the effect of genuine higher twists, we study the ratio of TMD to CGC near the back-to-back configuration k ⊥ ≈ 0, where kinematic power corrections vanish (see Sec 3.4). In Fig. 7 we show this ratio in nuclear DIS, and observe that if either the relative momentum P ⊥ or the virtuality Q are large enough this ratio goes to unity, and thus the genuine twists are suppressed. However, significant differences are observed at low to moderate values of Q and P ⊥ . They amount to more than a factor of 2 difference between the CGC and the TMD framework. In particular, a suppression of the CGC relative to the TMD prediction is observed at all back-to-back kinematics in the case of transverse polarization, while both suppression and enhancement are seen for the longitudinally polarized case depending on the virtuality and relative momenta.
The genuine saturation contribution to dijet production in the back-to-back limit k ⊥ → 0 (and in the light-cone gauge) is the result of multiple scattering, such that the total transverse momentum transferred to the qq pair is zero (since it compensates between several scatterings). The TMD and ITMD calculations miss this contribution as their result is proportional to the number of gluons with zero transverse momentum

Elliptic Anisotropy
In this section we study the elliptic anisotropy in the angle φ separately for transversely and longitudinally polarized photons. This quantity can be accessed by measuring the distribution of the momentum imbalance k ⊥ with respect to the relative momentum P ⊥ .
In the TMD framework, the magnitude of the elliptic anisotropy is proportional to the ratio of polarized to unpolarized Weizsäcker-Williams TMD (see Eqs. (3.60) and (3.61)): The study of v 2 as a function of the dijet momentum imbalance k ⊥ is a promising observable to extract the behavior of the linearly polarized gluons inside protons and nu-clei. In particular, it is expected that xh 0 /xG 0 → 0 in the regime of small momentum imbalance (k ⊥ Q s ), while the gluons are completely linearly polarized xh 0 /xG 0 → 1 for the perturbative limit k ⊥ Q s [72][73][74]. However, it is important to note that kinematic power corrections in the ITMD hard factor also induce correlations between the P ⊥ and k ⊥ (see Eqs (3.70) and (3.71)). These azimuthal correlations also couple to the linearly polarized WW. Consequently, distinguishing contributions purely from kinematic twists and those solely from the linearly polarized WW is not possible. However, we can compare the full ITMD predictions to the case where we set the linearly polarized WW gluon distribution to zero by hand, which will show the contribution from kinematic twists only. The comparison is shown in Fig. 8 as a function of the momentum imbalance k ⊥ . For both polarizations we find that the kinematic power corrections result in positive contributions to v 2 , which grow with the dijet momentum imbalance. On the other hand, in the TMD curve (which does not include kinematic power corrections) the elliptic anisotropy is sourced only by the linearly polarized gluons; in this case the sign of v 2 depends on the polarization of the photon (see Eqs. (3.60) and (3.61)). This additional v 2 from kinematic twists should be taken into account when extracting the linearly polarized gluon distribution in azimuthal dijet measurements [38]. In addition to kinematic corrections, genuine higher twists contribute to azimuthal correlations through higher body operators (see Sec. 3.4). Our goal in this section is to quantitatively study the different contributions to v 2 by comparing CGC, ITMD and TMD formalisms.
In Fig. 9 we present our results for the elliptic anisotropy as a function of the momentum imbalance k ⊥ when the virtual photon in DIS is transversely polarized. We consider proton and gold targets and different P ⊥ and Q 2 . We observe sizeable deviations between the TMD and CGC results when P ⊥ ∼ k ⊥ , signaling the appearance of azimuthal correlations due to kinematic power corrections and genuine higher twists. In particular, in the CGC we observe minima in v 2,T , which are absent in the TMD limit, and the anisotropies are considerably reduced at large k ⊥ (as argued previously the contribution to v 2,T from the intrinsic correlations due to linearly polarized gluons and those from kinematic twists come with different signs, see Fig.8).
In order to distinguish the relative importance of kinematic and genuine saturation corrections, we compare the CGC to the ITMD framework. We observe that the deviations in v 2T are enhanced for the nuclear DIS and for the lower scales in P ⊥ and Q. Interestingly, our results also indicate that the contributions of kinematic and genuine saturation contributions appear to come with different signs. In particular, we observe that at low momentum imbalance k ⊥ this results in an apparent better agreement between the CGC and its TMD limit despite the exclusion of kinematic power corrections.
Our results for the elliptic anisotropy for longitudinally polarized photons are shown in Fig. 10. In this case we observe that kinematic power corrections in the ITMD increase the azimuthal modulations compared to the TMD limit, which results in very large v 2,L at large momentum imbalance. In the CGC, this enhancement is partially tamed by genuine saturation corrections (most pronounced at low values of P ⊥ and Q 2 and for nuclear DIS).
Note that in the TMD framework v 2,L cannot exceed 1/2, which is because of the bound in Eq. (3.58). On the other hand, in ITMD and CGC the anisotropies exceed 1/2 due to kinematic (k ⊥ /Q ⊥ ) and genuine (Q s /Q ⊥ ) power corrections. Since the cross-section must remain positive, this indicates the presence of higher azimuthal correlations, which render the cross-section finite as we will show in the next section.

Quadrangular Anisotropy
In this section we study the quadrangular anisotropy v 4 in the azimuthal angle between P ⊥ and k ⊥ . We note that this contribution is completely absent in the TMD framework. Contributions to the quadrangular anisotropy have been computed for the first time in [83] by considering the first correction to the leading power of the dipole size r ⊥ when taking the TMD limit. Such higher powers of the dipole size can account for k ⊥ /Q ⊥ and Q s /Q ⊥ corrections without distinction unless one is careful with the separation between kinematic and genuine higher twists. This distinction for v 4 is the purpose of this section, where our results will be presented in the ITMD and CGC frameworks.
A quadrangular anisotropy can be generated by correlations between P ⊥ and k ⊥ in the hard factor in the ITMD (see Eqs. (3.70) and (3.71)) as well as correlations in the genuine higher twists embodied in Eq. (3.75). Fig. 11 demonstrates in scattering on Au targets that while the linearly polarized WW gluon distribution produces only an elliptic anisotropy in the TMD framework, in the ITMD framework the linearly polarized WW gluon distribution also has a contribution to v 4 , as it couples to higher modes. To illustrate this, observe that using Eqs. (3.57) and (3.68), we can write the differential cross-section for the ITMD as  Figure 11. Quadrangular anisotropy in γ * λ + Au → q +q + X: for transversely polarized photons (left), and longitudinally polarized photons (right). We show results in the TMD and ITMD framework. To illustrate the effect purely from kinematic twists, we also show the result for the ITMD in which we turn off the linearly polarized WW gluon TMD (dashed dotted).
where we introduced the hard factors: One could then expand H G and H h in modes in φ: Then we can obtain expressions for the elliptic and quadrangular anisotropies: .

(5.13)
This implies that when H h 4 (P ⊥ , k ⊥ ) is non-zero, xh 0 (x, k ⊥ ) can contribute to v 4 . While analytic expressions for these modes in the ITMD are difficult to obtain, the numerical results in Fig. 11 suggest that H h 4 (P ⊥ , k ⊥ ) is larger than H G 4 (P ⊥ , k ⊥ ). Turning off the linearly polarized WW gluon TMD by hand in the ITMD framework reveals the v 4 purely from kinematic twists. Its sign is positive for both transverse and longitudinal photon polarizations, while the full ITMD result has opposite signs for the two polarizations.
We present the CGC and ITMD results for the quadrangular anisotropy as a function of momentum imbalance k ⊥ at two different values of P ⊥ and Q 2 in Fig. 12 for the case in which the virtual photon is transversely polarized. First, we note that the quadrangular anisotropy v 4,T is negative (similar to v 2,T ) and on the order of a few percent at low momentum imbalance (k ⊥ P ⊥ ). In the lower panels (higher P ⊥ and Q) we observe that the results between the ITMD and CGC schemes closely agree with each other, showing that the quadrangular anisotropy is mostly driven by kinematic power corrections. The genuine higher twists tend to suppress v 4,T , this effect can be seen more clearly in the upper right plot corresponding to nuclear DIS and at low values of P ⊥ and Q.
Our results to the quadrangular anisotropy when the photon is longitudinally polarized are shown in Fig. 13. The behavior of v 4,L is similar to that of v 2,L : It is positive and significantly larger in magnitude than in the transversely polarized case. We observe that the quadrangular anisotropy is mostly sourced by kinematic power corrections; except for nuclear DIS and at small Q and P ⊥ (upper right panel). As in the transversely polarized case, the effect of the genuine higher twists is to suppress the quadrangular anisotropy.
We close this section by pointing out that our results for v 4,L and v 4,T as a function of momentum imbalance k ⊥ and at P ⊥ = 4 GeV and Q 2 = 10 GeV 2 have the same sign and similar magnitudes to those estimated in [74]. However, for this choice of kinematics the ITMD seems to be sufficient to describe the quadrangular anisotropy, especially in γ * + p scattering, and thus one only needs the WW gluon TMD (supplemented with ITMD hard factors) and no higher order derivative operators of Wilson lines.

Conclusions
The uncovering of the small-x structure in protons and nuclei is one of the major goals of the Electron-Ion Collider. The study of azimuthal dihadron and dijet anisotropies can shed light both on the emergence of gluon saturation as well as the nature of the Weizsäcker-Williams gluon TMD at small x (both linearly polarized and unpolarized). In this manuscript, we explicitly computed the production of quark anti-quark dijets at leading order in the ITMD scheme and the CGC EFT, extending the results in [1]. We briefly reviewed the origin of kinematic twists and genuine higher twist corrections by expressing the product of Wilson lines as a transverse gauge link at x − = −∞ [82], generalizing the small dipole size expansion in [70]. The genuine higher twists (genuine saturation effects) enter with powers of Q s /Q ⊥ , and these are present in addition to the saturation contributions Q s /k ⊥ in the WW gluon TMD at small x. Meanwhile, the ITMD resums kinematic power corrections to all orders in k ⊥ /Q ⊥ .
We found that for kinematics where either P ⊥ (dijet relative momentum) or Q (virtuality of the DIS photon) is sufficiently larger than the saturation scale Q s , the ITMD framework provides a good approximation to the CGC quark anti-quark differential crosssection and its anisotropies in the azimuthal angle φ between P ⊥ and k ⊥ . Compared to the standard TMD framework, we observe that the resummation of kinematic power corrections in the ITMD hard factors result in additional azimuthal correlations between P ⊥ and k ⊥ ; in particular, we observe a non-zero quadrangular anisotropy. Thus, we expect that a proper extraction of the linearly polarized WW gluon TMD xh 0 (x, k ⊥ ) from dijet azimuthal asymmetries [38] will require the implementation of the ITMD hard factors that have been computed analytically in this manuscript. An advantage of the ITMD framework is that it only depends on the WW gluon TMD in contrast to the CGC, which requires the computation of the quadrupole. This drastically simplifies the complexity of the computation of the differential cross-section for a wide set of kinematics, and makes it suitable for coupling to Monte-Carlo based hadronization and fragmentation routines. Efforts in this direction can be found in [123] where projections for the EIC have been made using the ITMD * 15 framework coupled to the event generator in [124].
On the other hand, genuine higher twists effects which are absent in the ITMD framework, are most significant for scattering off a heavy target, as expected due to the enhancement of the nuclear saturation scale Q 2 s ∼ A 1/3 . The genuine higher twist contributions suppress both the differential yield (near the back-to-back configuration) as well as the anisotropies (elliptic and quadrangular) in the angle φ for realistic EIC kinematics 16 . These genuine saturation effects may be accessed in the back-to-back measurement of low transverse momentum dihadrons and at low photon virtualities, and could be an important experimental tool for the gluon saturation searches at the EIC [37]. For non back-to-back configurations, it might be possible to probe genuine saturation effects when the virtual photon is longitudinally polarized as Fig. 6 showed deviations between CGC and ITMD across a wide range in k ⊥ . In Fig. 7 we observed that in e + Au collisions the back-to-back peak for quark anti-quark production in the CGC is suppressed by a factor of 2 relative to the TMD framework for P ⊥ , Q ∼ 1 GeV, whereas when either P ⊥ or Q ∼ 4 GeV both frameworks agree to a less than a few percent difference.
To go towards more phenomenological applications, we plan to incorporate in our analysis the effect of Sudakov resummation [125,126] and final state soft gluon radiation [127,128], which have been shown to significantly impact the measurement of azimuthal dihadron and dijet anisotropies at the EIC [37,123,129]. Furthermore, recent progress towards next-to-leading order computations for DIS at small-x [58, 105-107, 109-111, 130-133] will allow us to extend our computation to higher accuracy (for a recent computation of NLO contributions for dijet production in the CGC see [134]). The inclusion of parton showers, hadronization, full jet reconstruction, and appropriate background processes will be necessary to fully assess the impact of our study on future EIC phenomenology. Finally, it would be interesting to compare the results of the suppression of the back-to-back peak due to gluon saturation with those caused by the momentum broadening due to cold nuclear matter energy loss and coherent power corrections [135] (see also [136] for a recent computation in the CGC). 15 The ITMD * framework is equivalent to the ITMD framework in the absence of linearly polarized gluons, see Appendix E. 16 Except for the angle averaged differential yield for the longitudinally polarized photon, where we observe an enhancement of the CGC relative to the TMD when the dijet momenta are large and the virtualities are small (see Fig. 7).

A Wilson lines and transverse gauge links
The goal of this section is to briefly summarize the derivation of Eq. (2.3) which relates the product of two Wilson lines to a transverse gauge link at x − = −∞: As discussed in Sec. 2 this relation is helpful to establish the relation between CGC and TMD amplitudes rigorously since it allows for a systematic expansion in powers of gÃ ⊥ . More details and generalizations of the presented derivations can be found in [78,79,82].
To begin we recall that the gauge transformation (introduced in Eq. (3.14)) going from ∂ µ A µ = 0 Lorenz gauge toÃ + = 0 light-cone gauge is given by Similarly, the product of two Wilson lines can be expressed as the product of two gauge rotations at the x − = −∞ boundary In order to connect this result with that in Eq. (A.1) we shall express the product of these gauge rotations as a transverse gauge link. This follows by first noting where in the second equality we related the derivative of the gauge rotation to the gauge field inÃ (3.15)). From Eq. (A.5) then we find the recursive relation and in a similar fashion one can show These two recursive relations are satisfied by The integral is independent of the choice of path connecting the transverse points y ⊥ and x ⊥ due to the non-Abelian Stokes' theorem and the fact that the transverse components of the field strength are zeroF ij = 0 (at a given x − ,Ã i ⊥ (x − , x ⊥ ) is a pure gauge field in two dimensions), more precisely P exp C dz i ⊥Ã i ⊥ (z ⊥ ) = 1 for any closed path C contained in the transverse plane 17 .
It follows from Eq. (A.6) that the product of two Wilson lines satisfy an identical recursive relation After the application of the recursive relation in Eq. (A.7) we find where we define We would like to obtain a tractable expression for Eq. (A.11) which will be the starting point for our study of the production of quark and anti-quark pairs in the ITMD framework. where r ⊥ = x ⊥ − y ⊥ . Then we have The integral over s is easily done in momentum space using the identitỹ 15) and the simple integral over the phase The non-Abelian Stokes' theorem reads P exp C dxµA µ = P exp S dσµν U F µν U † where dσµν is the surface measure on S and U denotes a Wilson line connecting the point x ∈ S enclosed by the surface measure to an arbitrary base point O on C (see [137]).
We find Despite the fact that this expression might look more formidable than Eq. (A.11), we will show in Sec. 3.3 that the resultant differential cross-section for quark and anti-quark dijet production has a relatively simple analytic expression. Following the same procedure and with some algebra we find a similar expression for Eq. (A.12)

B Useful integrals
We list some useful transverse integrals: (B.6)

C Explicit representation of Dirac spinors
We work in the Dirac basis for gamma matrices where 1 is the two-by-two identity matrix.
The helicity operator h is defined as It is straightforward to check that the (massless) Dirac equation has the following solutions 18 where the subscripts ± denote the helicities 19 , φ k is the azimuthal angle of k ⊥ , and the normalization is chosen so that where the barred spinors are defined as usual byū = u † γ 0 . The following identities hold: D Computing the perturbative factors N λ To compute N λ σσ (r ⊥ ) we first perform the l − and l + integrations. The former is done immediately using the delta function δ(k − 1 − l − ), while the latter is performed via contour integration using Cauchy's theorem, we find To compute the transverse integration, we must first work out the Dirac algebra. We find where we used Eqs. (C.6) and (C.7), and ij λ,i ⊥ = −iλ λ,j ⊥ . This last identity of the polarization vector holds for circularly polarized states ±1 ⊥ = 1 18 In the massless case, the spinors corresponding to particle and anti-particle are the same, but correspond to opposite helicities. 19 Note that in this manuscript, p − is the large component of the spinor momenta; thus p 3 < 0. and The first steps to perform these integrals explicitly is to integrate the r i ⊥ factors by parts into derivatives w.r.t. p i ⊥ , and to use the integral form e i(k ⊥ ·r ⊥ ) − 1 i(k ⊥ · r ⊥ ) = Let us now address the α integral. Introducing (E.9) and using the standard integrals we have: With the help of the relation , (E.14) and the following additional relation for J ij T : Along with tedious but straightforward algebra, we find eventually: , (E. 16) and Therefore, we find the following expressions for the ITMD hard factors in the amplitude of quark anti-quark production: J λ=0,i σ,σ (P ⊥ , k ⊥ ) = 2(z 1 z 2 ) 3/2 Qδ σ,−σ , (E. 19) and J λ=±1,i σ,σ (P ⊥ , k ⊥ ) = (z 1 z 2 ) 1/2 [(z 2 − z 1 ) + σλ] δ σ,−σ (E.20) It is illustrative to consider two interesting limits: the so-called ITMD * limit [91], and photo-production limit Q 2 → 0.
The hard factors in the ITMD * scheme can be obtained in a diagrammatic approach with off-shell gluons [80], and they can be obtained from our results by the following projection [91]: Therefore, in the ITMD * limit it is enough to consider the projections: k j ⊥ J λ,j σσ (P ⊥ , k ⊥ ) = d 2 r ⊥ 2π e ik 2⊥ ·r ⊥ − e −ik 1⊥ ·r ⊥ N λ σσ (r ⊥ ) .

(E.24)
We refer the reader to [91] where the differences between the ITMD * and ITMD schemes have been numerically studied for the electro-production of heavy quarks. Finally, in the photo-production limit, we find which is compatible with the result from [78].

F Operator definition of the WW gluon TMD
In this appendix we briefly review the relation between the operator definition of the Weizsäcker-Williams distribution and the definition in Eq. (3.54), which was fleshed out for the first time in [70]. In what follows let us go back to the general case where we have not fixed the gauge condition for the gauge field A µ,a . The operator definition of a generic gluon TMD is constructed from the bilocal operator field strength tensors appropriately dressed by gauge links for gauge invariance. Let us define the finite path Wilson lines along the light-cone direction as and that along the transverse direction (where the path is a straight line) as The Weizsäcker-Williams distribution is defined as which involves the future pointing staple shaped infinite gauge link: To make the connection of the operator definition in Eq. (F.3) with that given in Eq. (3.54), let us now work in the Lorenz gauge ∂ µ A µ = 0 and in the small x limit, in which case A i ⊥ is suppressed. We note that the (infinite) light-like Wilson line defined in Eq. (2.1) can be expressed as 5) and the derivative ∂ i acting on the Wilson line is given by where we used A µ = δ µ+ A + and thus F i+ (b − , b ⊥ ) = ∂ i A + (b − , b ⊥ ) in light-cone gauge. Thus we find where we usedÃ i Then it follows from Eq. (3.54): (F.8) where we used [b ⊥ , b ⊥ ] +∞ = 1 in ∂ µ A µ = 0 gauge with the subgauge condition that transverse fields cancel at +∞, which is allowed for both components in the small x limit. The equivalence to the operator definition in Eq. (F.3) follows (strictly speaking) in the limit x → 0, and by noting the relation between the CGC average and the operator expectation value: