Back-to-back inclusive dijets in DIS at small $x$: Gluon Weizs\"acker-Williams distribution at NLO

In JHEP 11 (2022) 169, we performed the first complete computation of the back-to-back inclusive di-jet cross-section in Deeply Inelastic Scattering (DIS) at small $x_{\rm Bj}$ to next-to-leading order (NLO) in the Color Glass Condensate effective field theory (CGC EFT). We demonstrate here that for di-jets with relative transverse momentum $P_\perp$ and transverse momentum imbalance $q_\perp$, to leading power in $q_\perp/P_\perp$, the cross-section for longitudinally polarized photons can be fully factorized into the product of a perturbative impact factor and the nonperturbative Weizs\"{a}cker-Williams (WW) transverse momentum dependent (TMD) gluon distribution to NLO accuracy. The impact factor can further be expressed as the product of a universal soft factor which resums Sudakov double and single logs in $P_\perp/q_\perp$ and a coefficient function given by a remarkably compact analytic expression. We show that in the CGC EFT the WW TMD satisfies a kinematically constrained JIMWLK renormalization group evolution in rapidity. This factorization formula is valid to all orders in $Q_s/q_\perp$ for $q_\perp, Q_s \ll P_\perp$, where $Q_s$ is the semi-hard saturation scale that grows with rapidity.


Introduction
The Electron-Ion Collider (EIC) [2] will provide fundamental insight into the spatial and momentum distributions of partons in protons and nuclei. A particular promise of this high energy collider is the likelihood of robust extraction of transverse momentum dependent distributions (TMDs) and generalized parton distributions (GPDs) of gluons, with the latter being sensitive to spatial impact parameter distributions inside protons and nuclei [3,4]. At small x Bj , these distributions also provide a window into the onset of gluon saturation, characterized by an energy-dependent emergent saturation scale Q s , larger than intrinsically nonperturbative QCD scales [5,6]. A particularly intriguing prospect is that of extracting information on the gluon Weizsäcker-Williams (WW) TMD distribution at small x Bj . In light cone gauge, it represents the light cone quantized number distribution of gluons as a function of their transverse momenta in a proton or nucleus [7][8][9]. It was shown some time ago [10,11] that the leading order (LO) DIS back-to-back dijet cross-section, for dijets with relative transverse momentum P ⊥ and transverse momentum imbalance q ⊥ satisfying P ⊥ q ⊥ , can be expressed as the factorized product of a perturbative "impact factor" and the nonperturbative WW distribution. The WW distribution can be decomposed into two Lorentz structures: a trace component corresponding to the unpolarized WW distribution and a traceless component representing a linearly polarized distribution [12]. A first study of the behavior of the two distributions at LO in the CGC EFT [13][14][15][16][17][18] within the framework of the McLerran-Venugopalan model suggested that while these distributions behave similarly for q ⊥ > Q s , they have qualitatively different behaviour in the saturation regime of q ⊥ ≤ Q s [19,20]. A numerical study of azimuthal moments of the LO back-to-back distribution shows that the linearly polarized WW distribution is strongly suppressed for large nuclei in EIC kinematics for low q ⊥ but grows rapidly with increasing q ⊥ [21][22][23][24][25].
For robust predictions of high energy cross-sections in perturbation theory, it is essential that the computations be performed to at least next-to-leading order (NLO) accuracy that include corrections of O(α s ) to the LO cross-section . For back-to-back dijets, such NLO computations are essential to confirm that the aforementioned LO factorization into the impact factor and the WW TMD distribution holds to this accuracy. An additional essential feature of such NLO corrections is that they can be significantly enhanced in certain kinematic regions of the phase space of interest.
Of particular interest to us are "slow" NLO real and virtual gluon emissions, with light cone longitudinal momenta k + P + , where P + is the (large) momentum of the target nucleus. Their contributions are parametrically of magnitude α s ln(1/x) (x = k + /P + ) and are therefore O(1) contributions of the same magnitude as the LO contribution at small momentum fractions x 1. Further, there are contributions α N s ln N (1/x) ∼ 1 that appear at each N th loop order in perturbation theory. Fortunately, the all-order resummation of such leading logarithmic corrections in x (LLx) can be performed, ensuring the validity of perturbation theory which would otherwise break down. The renormalization group (RG) equation through which this is achieved is the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [52,53]. In the gluon saturation regime, due to the large power corrections from coherent multiple scattering that are enhanced by the saturation scale, the corresponding LLx equations generalize to the Balitsky-Kovchegov (BK) [54,55] and Balitsky-JIMWLK equations [54,[56][57][58][59]. The state-of-the-art are RG equations that resum next-to-leading-log α N +1 s ln N (1/x) contributions (NLLx); these are respectively the NLLx BFKL equation [60,61] and the NLLx JIMWLK equation [62][63][64].
These small x resummations are crucial to understand the energy (or rapidity) evolution of expectation values of correlators of Wilson lines that are the fundamental nonperturbative quantities in the CGC EFT; the WW distribution corresponds to one such correlator. In isolating such slow gluon contributions that first appear at NLO, one has to in practice introduce a cutoff x f to separate the large logs α s ln(x 0 /x f ) ∼ 1 (with x 0 denoting fast target modes k + ∼ P + ) from logarithms α s ln(x f /x Bj ) that are parametrically of order O(α s ) that are then absorbed into the impact factor. Indeed it is the required independence of cross-sections from this cutoff that leads to the small x RG equations which resum the small x logs to all orders in perturbation theory.
Since the small x RG equations are well-known, much of our focus in our previous papers [1,37] was on the computation of the NLO impact factor for dijets. For back-toback dijets, as discussed in [1], this computation of the impact factor is complicated by the appearance of so-called Sudakov logarithms that suppress the cross-section [65][66][67]. At NLO, these are the double logarithms α s ln 2 (P 2 ⊥ /q 2 ⊥ ) and single logarithms α s ln(P 2 ⊥ /q 2 ⊥ ) which can be of O(1) when the back-to-back dijet relative transverse momentum P ⊥ q ⊥ . The resummation of such large logarithms in perturbative QCD is achieved through the Collins-Soper-Sterman (CSS) formalism [68][69][70] (see also [71,72]). Numerical estimates of the impact of the resulting Sudakov suppression on semi-inclusive dihadron and dijet final states in small x kinematics accessible at the EIC suggest that such contributions are likely important and need to be quantified [73][74][75].
The first NLO study of the combined contributions from small x logs (from modes x 0 > x ≥ x f ) and the Sudakov logs from the impact factor (for x Bj < x ≤ x f ) was performed in [76,77] where the authors considered a number of final states in protonnucleus and electron-nucleus collisions, including back-to-back dijets. In particular, the factorization of double Sudakov logs from small x evolution was discussed in this work. Single Sudakov logs in this context were considered in [78][79][80]. In [1], we exploited the formalism for high order computations in DIS at small x developed in [32,33,81], and applied to inclusive dijets in [37], to compute all real and virtual contributions to inclusive back-to-back dijets. This yielded all Sudakov double and single logarithmic contributions (including terms not previously accounted for in the small x literature) as well as all finite O(α s ) contributions to the NLO impact factor computed for the first time.
A key finding of our work in [1] (also deduced simultaneously in [38] for photo-produced dijets) is that it is essential to impose lifetime ordering in addition to rapidity ordering on small x evolution in order to recover the right sign of the contribution from the Sudakov logarithms. Lifetime ordering was understood to be important previously but only in the context of modifying the kernel of the BFKL renormalization group (RG) equation [52,53] for fully inclusive cross-sections at next-to-leading logarithmic accuracy as observed for NLLx BFKL [60,61] in [82][83][84] and for the non-linear JIMWLK RG equation [56][57][58][59] in [85][86][87].
Significant insight is gained from the azimuthal angle moments cos(nφ) where φ is the angle between P ⊥ and q ⊥ of the back-to-back cross-section. At leading order, the zeroth moment of the cross-section is proportional to the unpolarized WW TMD distribution while the second moment is proportional to the linearly polarized TMD distribution. (All odd moments vanish due to quark-antiquark symmetry in the cross-section.) At NLO, due to their interplay with soft gluon radiation, both TMDs contribute to all even moments. In [1], we classified the NLO coefficient function -defined as all the NLO terms in the impact factor not containing Sudakov logarithms -as arising from two classes of diagrams. One class of such diagrams contain the contributions of virtual graphs in the back-to-back crosssection that manifestly satisfied TMD factorization for the correlation limit corresponding to P ⊥ q ⊥ . The other class of diagrams, corresponding to real and virtual diagrams where the gluon crosses the shockwave, appeared to have a complicated color structure different from the WW gluon TMD. When both the correlation limit and the dilute limit (P ⊥ q ⊥ Q s ) are satisfied, TMD factorization is manifestly recovered. However, when q ⊥ ∼ Q s , the color correlators did not appear to reduce to the WW TMDs in the correlation limit and therefore appeared to violate TMD factorization [1,38].
In this paper, we will revisit the computation in [1] of the NLO coefficient function contributing to the NLO impact factor 1 and demonstrate that TMD factorization can be proven to NLO when q ⊥ , Q s P ⊥ . The terms that previously appeared to violate TMD factorization in the dense saturated regime of q ⊥ ∼ Q s (real and virtual diagrams corresponding to gluons crossing the shockwave) are shown explicitly to satisfy TMD factorization. We compute the coefficient functions, including so-called "hard factors" and other finite O(α s ) terms, for all the factorized contributions. For the longitudinally polarized cross-section that we focus on in this paper, we obtain a remarkably compact and numerically tractable expressions for the coefficient functions in the factorized NLO result for angular moments cos(nφ) of the back-to-back cross-section.
We will present numerical results that demonstrate the relative magnitude of the contributions of the soft factor and those of the coefficient function on the NLO impact factor and their overall magnitude compared to the LO results. Detailed numerical studies including their implications for the extraction of WW distributions, and their interplay with the effects arising from gluon saturation, for both longitudinally and transversely polarized virtual photon-nucleus back-to-back dijet cross-sections, will be presented separately.
The paper is organized as follows. In Section 2, we give a brief summary of the main results of [1] and fix the notations and conventions used throughout the paper. The third and fourth sections are dedicated to the analytic calculation of the back-to-back limit of virtual and real NLO corrections. In Section 5, we put everything together to obtain the final result which is given in Section 5.2 for the TMD factorized form of the back-to-back NLO dijet cross-section, with analytic expressions provided for the coefficient function and the soft factor. We also present in Section 6 a preliminary numerical study of the NLO results that demonstrates the relative importance of the different contributions. We end in Section 7 with a summary and a discussion of future work. The paper is supplemented by four appendices. Appendix A provides details of the computation of the "master integrals" that appear in the NLO hard factor from the virtual diagrams where the gluon crosses the shockwave. Appendix B provides a mathematical justification of the correlation expansion that shows that higher order terms are suppressed in the dijet back-to-back limit. Appendix C is dedicated to the analytic computation of the Sudakov form factor using the running of the strong coupling at two loops. Finally, appendix D collects some useful integrals.
2 Review of computation in [1] of inclusive back-to-back dijet crosssection In this section, we will summarize the main results of [1] for the inclusive back-to-back dijet cross-section in DIS at small Bjorken x Bj . This section is structured as follows. After kinematic preliminaries, we will begin with the result for the cross-section in the CGC EFT at LO, which is factorized as the product of the hard LO impact factor and the Weizsäcker-Williams gluon transverse momentum distribution. We will then discuss our NLO results. We first outline all the real and virtual contributions that appear at NLO and state the results for the Sudakov form factor that are explicitly derived in [1]. A key result is that recovering the right Sudakov factor requires imposing a kinematic lifetime ordering constraint on small x evolution. In addition to the Sudakov terms, there are additional contributions whose structure we will outline. In particular, we show that the lifetime ordering constraint is also required for the NLO coefficient function. The explicit computation of these contributions will be the primary focus of this paper. For the DIS process γ * λ + A → qq + X, γ * λ is the virtual photon with polarization λ = L, T (longitudinal or transverse) and q andq label the two produced jets (or quarkantiquark pair at the parton level). The four momenta of the two produced jets are labeled k µ 1 and k µ 2 , and q µ is the four-momentum of the virtual photon. We will work in the dipole frame where the virtual photon propagates along the minus light cone direction while the proton or nucleus A, with 4-momentum P µ A propagates along the plus light-cone direction: Table 1 summarizes all the relevant kinematic variables for this process. Throughout the paper, we will work in the high-energy limit W 2 Q 2 , P 2 ⊥ . The full DIS cross-section is determined by convoluting the cross-section for γ * + A → qq + X with the longitudinal and transverse photon fluxes defined as where s is the center of mass energy squared per nucleon in the e + A collision, y = Q 2 /(s x Bj ) denotes the inelasticity (in the limit s m 2 n where m n is the nucleon mass), virtuality squared of the exchanged photon s = (P n + k e ) 2 nucleon-electron center of momentum energy squared W 2 = (P n + q) 2 nucleon-virtual photon center of momentum energy squared and α em is the electromagnetic fine structure constant. For fixed s, the e + A → e + qq + X cross-section is given by This expression is valid to all orders in the strong coupling constant. It is therefore sufficient to consider only the perturbative expansion in α s of the γ * + A → qq + X cross-section. In the NLO discussions in this paper, we will focus on the longitudinally polarized γ * L dijet cross-section alone, the expressions for which are more compact and easier to derive than that of the transversely polarized γ * T dijet cross-section. The results for the transversely polarized photon case, in combination with the results here, will be summarized in a forthcoming paper.

Dijet back-to-back limit at LO
After a brief overview of the inclusive dijet cross-section in the CGC EFT at leading order, we define the back-to-back kinematics and provide the TMD factorized expression of the CGC formula in this limit. We refer the reader to Section 2 in [1] for details.

Leading order cross-section in the CGC EFT
As discussed at length in previous work [32,37,81,88], weak coupling computations in the CGC EFT can be performed using standard Feynman diagrams techniques in the light cone gauge A − = 0, supplemented by dressed vertices corresponding to the coherent multiple scattering of a quark or a gluon inside the classical (over occupied) gauge fields of the target. These vertices are denoted, respectively, by cross and bullet symbols in Fig. 1 and Fig. 2. Figure 1. Leading order contribution to the amplitude for dijet production. The cross symbol on the quark and antiquark legs refers to the CGC quark effective vertex. The four-momenta of the quark and antiquark read k For the leading order production of a quark-antiquark pair in DIS at small x Bj , the only relevant Feynman diagram is the one shown in Fig. 1. After evaluating and squaring the Feynman amplitude associated with diagram 1, the fully differential leading order crosssection for inclusive production of two jets in the CGC can be written in the compact form 2 In this expression, e 2 f is the sum of the squares of the light quark fractional charges, and δ is an overall longitudinal momentum conserving delta function. The differential measure d 8 X ⊥ comes from the transverse coordinate integration contained in the CGC vertices and reads with x ⊥ (y ⊥ ) the transverse coordinate at which the quark (antiquark) crosses the shockwave in the amplitude (and similarly with prime coordinates for the complex conjugate amplitude). Differences of generic transverse coordinates a ⊥ and b ⊥ will be labeled r ab with The integrand in Eq. (2.5) is expressed as the convolution of a perturbative "impact factor" R λ LO corresponding to the QED splitting of the virtual photon into the quarkantiquark pair with the color correlator Ξ LO describing the interaction of the pair with the small x gluons of the target. At LO, the explicit expressions for these impact factors are respectively, for longitudinally and transversely polarized virtual photons. The effective virtualityQ is defined asQ 2 = z 1 z 2 Q 2 and K n (x) are the modified Bessel functions of the second kind and order n.
The color correlator Ξ LO is a CGC average of the product of Wilson lines at the projectile rapidity scale Y . Indeed, any observable O computed in the CGC EFT must be averaged over the large x color charge configurations in the target nucleus as where the color charge configuration ρ A is drawn from a stochastic gauge invariant weight functional W Y [ρ A ] with Y = ln(z), for a given typical fraction z of the projectile q − momentum corresponding to the kinematics of the observable O.
For inclusive dijet production in DIS, this correlator is given by where V (x ⊥ ), in the fundamental representation of su(3) whose generators t a are the Gell-Mann matrices, is the lightlike path ordered Wilson line This Wilson line factor resums to all orders the coherent multiple scatterings between the quark and the small x gluons in the target. The gauge invariant observables corresponding to these Wilson lines factors are the aforementioned expectation values of the dipole D and quadrupole Q operators: (2.14)

Back-to-back kinematics
The back-to-back limit of the LO inclusive DIS dijet cross-section is defined by the hierarchy of scales, where P ⊥ and q ⊥ denote respectively the relative transverse momentum of the dijet system and the dijet transverse momentum imbalance (see Table 1). The tree-level kinematics for the process γ * + g → qq, where g denotes a gluon in the nucleus wave-function, imposes the following constraint for the plus longitudinal momentum fraction of the gluon relative to P + n : where the second equality is valid in the Regge W 2 Q 2 . We have also introduced the dijet invariant mass M qq defined as At small x, as noted previously, the saturation scale Q s (x) which grows with decreasing x emerges as a semi-hard scale that plays an important role in the many-body dynamics. A priori, there is no constraint on the magnitude of Q s relative to q ⊥ or P ⊥ . In this work, we will assume that P ⊥ Q s , (2.18) without any assumption on the magnitude of the ratio Q s /q ⊥ ; saturation effects are large 3 when q ⊥ Q s . When q ⊥ Q s , back-to-back dijets probe the dilute or universal TMD regime which is insensitive to gluon saturation.

LO TMD factorized expression
The overlap between the CGC and TMD formalism for back-to-back inclusive dijets was discussed previously in [10,90]. This overlap becomes manifest when one performs an expansion in both q ⊥ /P ⊥ and Q s /P ⊥ and keeps the leading terms in this double expansion. Formally, one initiates this power counting with the change of variables, 20) in the leading order CGC cross-section: followed by an expansion of the color correlator Ξ LO for u ⊥ b ⊥ and u ⊥ b ⊥ : (2.22) After the necessary Fourier transform is performed, each additional power of u ⊥ or u ⊥ brings an extra power of 1/P ⊥ times a transverse scale of order q ⊥ or Q s . The latter appears via the additional transverse coordinate derivative of a Wilson line in the Taylor expansion of Ξ LO . Therefore the truncation in Eq.(2.22) of the expansion of the color correlator captures the leading term in the limit P ⊥ q ⊥ , Q s and the order of magnitude of the neglected terms is O(max(q ⊥ , Q s )/P ⊥ ).
The operator in the leading term in this power counting is nothing other than the small x Weizsäcker-Williams distributionĜ ij whose Fourier transform is the WW TMD distribution: One can rewrite the LO inclusive dijet DIS cross-section in back-to-back kinematics as the factorized expression 4 (2.25) Here the LO coefficient function in this case can be expressed in terms of a "hard factor" H λ,ij LO , defined to be A straightforward calculation using Bessel identities gives for the leading order hard factors, respectively for the longitudinally and transversely polarized virtual photon cross-sections.
A further simplification occurs when one decomposes this distribution further into a trace part (the unpolarized WW gluon TMD G 0 Y (q ⊥ )) and a traceless piece (the linearly polarized WW TMD h 0 Y (q ⊥ )) as Substituting the above expressions back into Eq. (2.21) then gives the very simple factorized expressions for the LO longitudinally and transversely polarized back-to-back dijet crosssections, respectively: By inspection, the azimuthal angle average of Eqs. (2.30) and (2.31) is only proportional to the unpolarized WW TMD. Conversely, the cos(2φ) moment of these cross-sections (the elliptic anisotropy) is proportional to the ratio of linearly polarized to the unpolarized distribution. As we will discuss, soft gluon radiation at NLO quantitatively modifies these LO results. For future reference, the LO hard factor will appear in the NLO cross-section in the form of the trace Due to the large gluon occupancy at small x, the WW gluon distribution G ij Y (q ⊥ ) is parametrically of order O(1/α s ) [13]. Thus the LO cross-section is parametrically O(1), and the NLO corrections are of order O(α s ).
2.2 Back-to-back limit at NLO: Sudakov suppression and kinematically constrained small x evolution The NLO inclusive dijet cross-section at NLO in general kinematics was first computed in the CGC EFT in [37]. All the Feynman diagrams that contribute at this order are shown in Fig. 2 At NLO, there is an additional CGC vertex (represented by the bullet symbol) associated with the coherent multiple scattering scattering of the radiated gluon off the shockwave [81,[91][92][93]. In order to obtain an infrared and collinear finite result, one has to introduce jets. The jet definitions that we consider throughout this paper are generalized k t algorithms, which all give the same NLO impact factor up to powers of R 2 corrections, with R the jet radius [94]. Figure 2. Feynman diagrams that contribute to the inclusive dijet cross-section at NLO. Top: Real gluon emission diagrams. Middle: Self energy diagrams. Bottom: Vertex correction diagrams. Diagrams obtained from q ↔q interchange are not shown, and will be labeled with an additional prime index. (For example, R 2 → R 2 .) Only diagrams in which the gluon does not scatter off the shockwave (represented by the red band) contribute to the Sudakov double and single logarithms.
In [1], we computed the back-to-back limit of our NLO results in [37] and expressed the results in terms of the moments of the azimuthal angle φ between between P ⊥ and q ⊥ . For the zeroth moment of the distribution (the angle averaged cross-section), we obtained the result (up to terms of order O(R 2 )), For the second moment of the distribution 5 , we obtained Results for the higher moments, and for transversely polarized photons were also provided in [1] in the appendices. Final results for transversely polarized photons follow analogously to those that we will obtain here; as noted earlier, the combined results will be discussed in a subsequent paper. We will show in this paper that the expressions provided above in Eqs. (2.33) and (2.34) can be greatly simplified enabling one to numerical evaluate with relative ease. However before we proceed to do so in the following sections, it will be useful to first explain key features of the above results.
Sudakov form factor. The function S(P 2 ⊥ , r 2 bb ) is the Sudakov or "soft" form factor. At LO in the coupling, this quantity is unity and one recovers in the above expressions the result for the zeroth and second azimuthal moment of Eq. (2.30), respectively, in terms of the unpolarized and linearly polarized WW TMD distributions.
At NLO, the function contains the Sudakov double and single logarithms derived explicitly in [1], which give O(1) contributions when P ⊥ q ⊥ . Assuming CSS evolution holds at small x, the all order contribution of the resummation of these Sudakov logarithms is given by The coefficients of the single Sudakov log that we extracted in [1] are given by up to powers of R 2 corrections. Here ∆η 12 ≡ η 1 − η 2 , and z f , Q f are factorization scales associated with the rapidity evolution of the WW gluon TMD discussed below.
Kinematically constrained rapidity evolution. An important result from the study in [1] is that consistent resummation of Sudakov logarithms to all orders requires that one employs a kinematically improved small x renormalization group (RG) evolution of the WW gluon TMD. This RG equation reads, in integral form, respectively, the initial rapidity of the evolution and the rapidity factorization scale. The transverse momentum scale Q 2 f should be taken of order max(P 2 ⊥ , Q 2 ) but was chosen to be arbitrary in [1]. We will express Q f in this paper in terms of the z f and the kinematics of the dijet system.
The kernel of the rapidity evolution of the WW gluon TMD (which we shall denote henceforth as K DMMX ) was computed explicitly from the Balitsky-JIMWLK hierarchy in [95]. The notation "⊗" is meant to remind the reader of the specific application of the universal JIMWLK kernel to the LO WW correlator, which generates additional correlators which are not WW-like. In other words, as shown in [95], the rapidity evolution of the WW TMD does not have a closed form in the dense saturation regime. The main difference of our result with respect to the evolution found in [95] lies in the Θ-function term which enforces lifetime (or k + g ) ordering of successive gluon emissions from the projectile, on top of the usual projectile momentum k − g ordering. Hence in the computation of Eqs. (2.33) and (2.34), the WW gluon TMDs should be evolved employing this kinematically improved RG evolution equation.
NLO coefficient function. The virtual diagrams SE 2 , V 2 and V 3 (plus their q ↔q counterparts) -in which the gluon loop does not interact with the shockwave -generate NLO hard factors that are distinct from H λ,ij LO and are denoted H NLO,1 and H NLO,2 in the formula Eq. (2.33). Exact expressions for these are given in the next section; their explicit computation is a key result of this paper. Other finite terms of O(α s ) (not containing Sudakov or rapidity logs) are also present and given explicitly in Eqs. (2.33) and (2.34). They come from unresolved real gluon emission inside jet cones or soft gluon radiation outside the jet cones and are discussed further in Section 4.1. In particular, the latter induces a sensitivity of the zeroth (second) moment to the linearly polarized (unpolarized) WW TMD which is a pure NLO effect, as shown in the fifth line of Eq. (2.33) (Eq. (2.34)). All of these terms, namely, NLO hard factors and finite terms from real gluon emissions, together comprise the NLO coefficient function.
The "other" contribution. The last term in Eq. (2.33) labeled is the contribution to the back-to-back dijet production of all the diagrams in which the emitted gluon interacts with the shockwave. These terms are particularly difficult to compute because they depend on a number of other color correlators besides the WW correlator. This is because a naive correlation expansion similar to Eq. (2.22) does not yield the WW gluon TMD. However, we will see that there exists an expansion for these contributions which correctly accounts for the leading power back-to-back limit and where the leading term is the WW gluon TMD in coordinate space. The explicit calculation of this "other" contribution is another key result of this paper. This is because it establishes the highly nontrivial result that the NLO back-to-back dijet cross-section, for P ⊥ q ⊥ , Q s can be fully factorized into a NLO impact factor and the WW gluon TMD, providing a complete NLO proof of TMD factorization for this process. More precisely, the "other" contribution also introduce new hard factors, labeled H NLO, 3 and H NLO,4 in the next section, that contribute to the NLO coefficient function.

NLO hard factors from virtual diagrams
This section is dedicated to the calculation of the hard factors H NLO,1 , H NLO,2 and those that we will see emerge from the virtual diagrams SE 1 and V 1 included in the dσ (0),λ=L other term 6 of Eq. (2.33). We will obtain fully analytic expressions for these hard factors in terms of logarithms and dilogarithms of the ratio Q/M qq . The analytic results simplify considerably the otherwise challenging numerical computation of the small x back-to back dijet cross-section at NLO.
On a more conceptual level, the calculation of the back-to-back limit of diagrams SE 1 and V 1 demonstrates that the full virtual cross-section can be factorized into a perturbative hard factor and the WW gluon TMD in back-to-back kinematics, which in combination with the results of the next section, shows factorization of the full NLO cross-section. As previously noted, this is a highly nontrivial result. The hard factors for the virtual diagrams in which the gluon scatters off the shockwave are meaningful if and only if the rapidity divergence that is subtracted is kinematically constrained. Thus this section provides supplementary evidence (to that provided by the "right sign" of the Sudakov double logs [1,38]) of the necessity to go beyond naive small x evolution and to use the kinematic constrained evolution discussed in the previous section. Further, our computation here enables us to fix the arbitrary scale Q f in our kinematically improved rapidity evolution equation for the WW TMD, and therefore complete our calculation of the Sudakov single log coefficients in [1]. 6 A similar analysis can be carried out for the term dσ (2),λ=L other in Eq. (2.34). We will not work this term out explicitly here but only quote the final result along with that for Eq. (2.33).

Final state vertex correction
We begin with the hard factor H NLO,2 arising from the final state vertex diagram V 3 in Fig. 2. The leading term in the correlation expansion of the CGC color correlator in this diagram gives the WW gluon TMD. In the back-to-back limit, the V 3 × LO (and its complex conjugate) contribution to the cross-section, with its rapidity divergence subtracted, factorizes as [1]: One notes that its contribution is suppressed at large N c . Thus this diagram does not contribute in the large N c approximation. The hard factor H λ=L,ij NLO,2 (P ⊥ ) associated with this diagram reads [1] In [37], we were unable to find a closed form analytic expression for the diagram V 3 and our final result was expressed in terms of the function J defined by (3.5) The main result of this section will be the explicit analytic evaluation of H λ=L,ij NLO,2 (P ⊥ ). As discussed in the Appendix J of [37], the −i prescription in the denominator, coming from the Feynman propagator of the quark or antiquark prior to the emission of the virtual gluon, is important to specify the integration contour since the denominator has a pole on the positive real axis of In the formula for H λ=L,ij NLO,2 , the last term proportional to is the subtraction term of the leading z g → 0 divergence of the integral, which ensures that the integrand inside the curly bracket is finite. We postpone the calculation of this term to the end of this subsection since its form can be inferred from its defining property of regularizing the slow (z g → 0) divergence. Let us move now to the more complicated pieces. First, given the form of R λ=L LO (u ⊥ , u ⊥ ), one can obviously factorize the u ⊥ integral from the u ⊥ one. The latter gives (cf Eq. (D.2)): (3.7) In the u ⊥ integral, the term proportional to K 0 (−i∆ V,3 u ⊥ ) does not contribute because the phase proportional to P ⊥ · u ⊥ cancels and one is left with an integral which only depends on u i ⊥ times a function of u 2 ⊥ . We are thus left with the computation of the integral where the (z g → 0) term is the slow subtraction term given in Eq. (3.6). As outlined in the beginning of this section, the integral I i 1 can be computed analytically, even though J does not have a closed analytic expression. This calculation is detailed in Appendix A.2. In terms of ∆ V3 ,Q V3 and K ⊥ defined by Eqs. (3.10) The imaginary part of the integral, which comes from the −i prescription in the contour integration for J , does not contribute to the total cross-section since the complex conjugate term is added to the hard factor in Eq. (3.2). In principle, the "c.c." term should be added at the cross-section level and not in the hard factor. It is alright to do so since the Fourier transform of the WW gluon TMD is real. If the Fourier transform of the WW gluon TMD were instead complex, this imaginary part of I i 1 would have contributed to the total cross-section.
Subtracting the slow gluon divergence. The final step in the calculation of H λ=L,ij NLO,2 is to perform the remaining z g integral over the virtual gluon longitudinal momentum. As explained above, this integral is convergent provided the z g → 0 logarithmic divergence is subtracted off. In principle, one should then compute the subtraction term given by Even though the u ⊥ integral can be performed, one can actually infer the final result from the z g → 0 limit of Eq. (3.10). Inserting the values of K ⊥ , ∆ V3 andQ V3 , we find that Hence, since by definition H λ=L,ij NLO,2,slow is equal to the leading slow divergence of H λ=L,ij NLO,2 , we have With this subtraction term, the hard factor can be expressed as a single integral over z g : with the variable χ defined as The integral over z g is now perfectly regular. After a rather tedious calculation, the z g integral can be expressed in terms of logarithms and dilogarithms as where the dilogarithm or Spence's function is defined by for x ≤ 1. Note that Eq. (3.16) is finite when z 2 − z 1 χ 2 = 0 thanks to the cancellation of the single pole (which relies on the identity z 1 + z 2 = 1), and that the arguments in the Spence's function are always smaller than 1.

Initial state self-energy and vertex corrections
The virtual diagrams SE 2 and V 2 also contribute to the back-to-back limit through another hard factor labeled H λ=L,ij NLO,1 (P ⊥ ) in Eq. (2.33), whose definition is [1]: To compute this hard factor, one notices that it is related to H λ=L,ij NLO,2,slow (P ⊥ ) since one has Using the result Eq. (3.13), one gets the simple expression where we recall that χ = Q/M qq . We have obtained this expression based on our knowledge of the slow gluon divergence of H λ=L,ij NLO,2 (P ⊥ ) but we note that one can also obtain it by a direct computation of the integral, as shown in Appendix A.1.

Self-energy crossing the shockwave
We turn now to the discussion of the back-to-back limit of virtual diagrams in which the gluon scatters off the shockwave, namely the self-energy SE 1 and the vertex correction V 1 diagrams. Their respective contributions are included in the term dσ (0),λ=L other because they depend on color correlators which do not "naturally" reduce to the WW gluon TMD. However, we will see upon closer inspection that their back-to-back limit does depend only on the WW gluon TMD, again up to q ⊥ /P ⊥ power corrections.
We first consider diagram SE 1 , the expression for which was obtained in [37]. For a longitudinally polarized photon, the diagram SE 1 times the complex conjugate of the leading order diagram reads where the color correlator Ξ NLO,1 is given by Note that this contribution to the cross-section depends on the integral over the transverse coordinate z ⊥ representing the position where the gluon crosses the shockwave. The formula for SE 1 also involves the qqg effective dipole size X V defined by The square bracket in Eq. (3.21) contains two terms whose difference regulates the UV divergence r zx → 0 of the self-energy. This counterterm is discussed at length in [37]. Finally, the slow divergence subtraction term regularizes the z g → 0 divergence of the z g integral and will be discussed towards the end of this subsection. In order to obtain the back-to-back limit of the expression Eq. (3.21), one must absorb the phase into an overall phase factor that resembles the leading order one with a suitable change of variables. One possible such change of variables is There is a natural physical interpretation for this change of variables in terms of the transverse location of the parton in the diagram. The transverse coordinate u ⊥ is the difference between the transverse locations of the quark after reabsorbing the gluon (and also before emitting it) x ⊥ = (1 − z g /z 1 )x ⊥ + z g /z 1 z ⊥ and the transverse location of the antiquark at the position y ⊥ where it interacts with the shockwave. In a similar fashion, the transverse location of the photon b ⊥ can be expressed as b ⊥ = z 1 x ⊥ + z 2 y ⊥ as in the leading order case Eqs. (2.19). Finally, the variable r ⊥ corresponds, up to the factor 1 − z g /z 1 , to the transverse separation between the quark and the gluon at the shockwave (see Fig. 3).
where Ξ NLO,1 represents the color correlator Ξ NLO,1 expressed as a function of u ⊥ , b ⊥ , r ⊥ , u ⊥ and b ⊥ . The variable ω appearing in one of the Bessel functions is defined to be Notice that in the UV subtraction term, we have used the LO change of variables given by Eqs. (2.19) as this term is not multiplied by the phase Eq. (3.24).
Correlation expansion of Ξ NLO,1 . In Eq. (3.27), it is not yet manifest that the back-toback limit P ⊥ q ⊥ of diagram SE 1 only depends on the WW gluon TMD since Eq. (3.27) depends on the color correlator Ξ NLO,1 . To see this, one must evaluate Ξ NLO,1 as a function of u ⊥ , b ⊥ and r ⊥ using Expanding Ξ NLO,1 , to get the leading power in q ⊥ /P ⊥ or Q s /P ⊥ , it is sufficient to expand the correlator to leading order in u ⊥ , u ⊥ and r ⊥ . In Appendix B we show that any higher order terms in the Taylor series are power suppressed in the ratio q ⊥ /P ⊥ or Q s /P ⊥ . Hence even though one integrates over r ⊥ , the back-to-back limit of the diagram SE 1 is controlled by values of r ⊥ such that r ⊥ ∼ u ⊥ ∼ 1/P ⊥ -namely, gluons with large transverse momentum. The Taylor series expansion of the color correlator Ξ NLO,1 around b ⊥ (dubbed the "correlation expansion" throughout this paper) gives to first finite order, One recognizes the WW color correlator in the second line. Hence as conjectured, the product of the diagram SE 1 with the LO diagram can be factorized into a hard factor and the WW gluon TMD: with the hard factor NLO, 3 defined as The term proportional to r i ⊥ u j ⊥ that appears in Eq. (3.30) does not contribute since it gives an integral of the form that vanishes for any function f . The rapidity subtraction term labeled H λ=L,ij NLO,3,slow (P ⊥ ) ensures that the z g integral is integrable near z g = 0 and is obtained by the same expansion of the color correlator in the rapidity subtraction term dσ λ=L SE 1 ,slow . As we shall see, this subtraction term is given by the kinematically constrained slow gluon limit of diagram SE 1 .

Calculation of H λ=L,ij
NLO,3 (P ⊥ ). Similarly to the cases of the NLO, 1 and NLO, 2 hard factors, the hard factor associated with SE 1 can be computed analytically. Since the correlation expansion assumes u ⊥ ∼ r ⊥ , one cannot make any approximation in the transverse coordinate integral over r ⊥ so we have to evaluate it as is. The u ⊥ integral is straightforward from Eq. (D.2): where the scalar integral I 2 is defined by Details of this result are given in Appendix A.3. We quote here the final result for H λ=L,ij NLO,3 before integrating over z g , which reads with χ = Q/M qq , as previously.
Subtracting the slow gluon divergence. The z g integral is double logarithmically divergent as z g → 0, but this divergence should be cured by the slow subtraction term dσ λ=L SE 1 ,slow in Eq. (3.27) that we have not specified so far. A first attempt would be to use the subtraction term used in [37] that enables one to obtain the first step in the JIMWLK evolution of the LO cross-section between z 0 and z f . This naive subtraction term is given by setting z g = 0 everywhere but in the 1/z g prefactor. One then gets: Expanding Ξ NLO,1 as above, we get for the hard factor This naive subtraction term is problematic for two reasons. Firstly, it results in the integral which is IR divergent when r ⊥ goes to infinity. It means that one is not allowed to expand the color correlator Ξ NLO,1 inside the r ⊥ integral. Secondly, even if we assume that the r ⊥ integral is naturally cut by a scale of order say 1/Q s , the subtraction term would be still proportional to a single logarithm of the z 0 cut-off, which prevents the complete cancellation of the double logarithmic divergence in Eq. (3.37). The solution to this problem was already outlined in Section 2.2 and points independently to the kinematically constrained rapidity evolution of the WW gluon TMD discussed there. Thus not only is the kinematic constraint crucial to obtain the correct sign for the double logarithmic Sudakov suppression of the back-to-back dijet cross-section as we showed previously in [1], it is also mandatory to obtain a well defined hard factor for the virtual graphs where the gluon scatters off the shockwave. This is a nontrivial consistency check of our results for back-to-back dijets in the CGC, but hardly a surprise. This is because the original motivation of the existence of a kinematic constraint followed from examining the spacetime structure of the diagrams SE 1 and V 1 as z g → 0 [1]. (See [96] for similar arguments in the context of fully inclusive DIS.) Indeed, for the NLO contributions in which the gluon crosses the shockwave, z g ∝ ω and r ⊥ are coupled in the NLO light-cone wavefunction: Thus, in order to recover the expected LO wavefunction K 0 (Qu ⊥ ) in the slow gluon limit z g → 0, we must require the auxiliary condition: which we refer to as the kinematic constraint. In more concrete and mathematical terms, the constraint Eq. (2.39) applied to the rapidity evolution of the WW gluon TMD translates into the improved subtraction term for the diagram SE 1 : where Q f is a priori an arbitrary transverse momentum scale of order P ⊥ , based on the argument leading to Eq. (3.43). Expanding Ξ NLO,1 to first non trivial order to extract its component proportional to the WW gluon TMD, we get the subtraction term of the hard factor H λ=L,ij NLO, 3 : (3.45) The integral over r ⊥ is now convergent both in the UV and in the infrared: The resulting integral over u ⊥ of the second term in this expression is given by Eq. (D.4) in the appendix. At the end of the day, we get We observe that the k − g or z g dependence of the kinematic constraint turns the single logarithmic slow divergence into a double logarithmic one. Therefore, the use of a kinematically constrained rapidity evolution solves both of the problematic issues mentioned after Eq. (3.40) resulting from "naive" leading logarithmic BFKL evolution. This is by no means the end of the story and indeed results in a further satisfying outcome. The counterterm H λ=L,ij NLO,3,slow (P ⊥ ) must cancel both double and single logarithmic divergences of H λ=L,ij NLO,3 (P ⊥ ) as z g goes to 0. This requirement for a convergent z g integral entirely specifies the scale Q f which was treated as an arbitrary scale in our previous work. Comparing Eq. (3.47) with Eq. (3.37) for H λ=L,ij NLO,3 (P ⊥ ), one finds that cancellation of the z g divergences requires that (3.48) Since typically z f ∼ z 1 z 2 , one sees that indeed, Q f must be of order max(Q, P ⊥ ). Of course, this identity is specific to our implementation of the kinematically constrained rapidity evolution of the WW gluon distribution in coordinate space. In this sense, it is a scheme-dependent statement. However the ratio z f /Q 2 f is scheme-independent since it only depends on the kinematics of the dijet process. Therefore the constraint Eq. (3.48) may lead to scheme-independent results, one of them being the value of the Sudakov single logarithmic coefficient s f defined in Section 2.2. We shall come back to this point in Section 5.1.
Final result. With the counterterm in Eq. (3.47), and the constraint Eq. (3.48), the third NLO hard factor Eq. (3.37) is now finite since the z g integral becomes regular for z g = 0. The calculation of this integral is straightforward and gives The integral in the second line is the leftover contribution from the subtraction of the slow gluon phase space between z f and z 1 . It does not need to be computed as it will cancel with the C F dependent slow gluon counterterm of diagram V 1 , as we shall see in the next section. Without this important cancellation, the integral in the second line would generate double logarithms of the rapidity factorization scale z f , spoiling rapidity factorization. Finally, the diagram SE 1 obtained by quark-antiquark exchange (z 1 ↔ z 2 ), and the complex conjugate terms must be added together to obtain the final self-energy contributions of gluons crossing the shockwave to the back-to-back dijet cross-section.
Order of magnitude of neglected terms. To summarize our discussion thus far, we showed that the leading term in the correlation expansion of SE 1 × LO can be factorized into a perturbative hard factor and a nonperturbative WW TMD distribution. Our method relied on (i) finding suitable transverse coordinate integration variables such that u ⊥ − u ⊥ is conjugate to P ⊥ and likewise b ⊥ − b ⊥ is conjugate to q ⊥ , as is the case for the LO crosssection, resulting in the coordinate phase space factor (ii) Expanding the CGC operators inside the transverse coordinate integration, to first order around x ⊥ = y ⊥ = z ⊥ = b ⊥ in the limit |r ⊥ |, |u ⊥ | |b ⊥ | where r ⊥ is typically conjugate to the transverse momentum of the virtual gluon. While it is natural to justify the small |u ⊥ | expansion from the phases Eq. (3.50), one may wonder why one is allowed to truncate the expansion to first order in r ⊥ as well, to arrive at the back-to-back limit. Hence to prove TMD factorization of all virtual diagrams in the leading power approximation, we must also demonstrate that the higher order terms in the correlation expansion of the CGC operator Ξ NLO,1 , are power suppressed either by q ⊥ /P ⊥ or Q s /P ⊥ .
As a matter of fact, this can be understood from a simple dimensional analysis argument. The kinematic constraint from the NLO light cone wave-function (the Bessel function) imposes r ⊥ and u ⊥ to be of comparable sizes, r ⊥ ∼ u ⊥ . Hence, after Fourier transform, each additional r ⊥ power leads to a 1/P ⊥ ∼ 1/Q suppression. We provide in Appendix B.1 a full mathematically rigorous derivation showing that this is indeed the case.

Vertex correction crossing the shockwave
The back-to-back limit of the vertex correction with gluon interacting with the shockwave, labeled V 1 in Fig. 2 can be addressed in a similar fashion. In particular, this diagram can also be factorized into a hard part and the WW gluon TMD. The most important part of this calculation will be to cross-check that the same value of Q f as the one extracted from diagram SE 1 in Eq. (3.48) must be used in the kinematically constrained rapidity subtraction term to ensure a finite hard factor for the vertex diagram V 1 .
This diagram was computed in full generality within the CGC EFT in [1]. For a longitudinally polarized virtual photon, we found for the product between V 1 and the complex conjugate of the leading order amplitude. As the phases that appear in this expression are identical to those present in Eq. (3.21) for SE 1 , we shall use the same change of variables given in Eq. (3.26). Another motivation for using the same values of u ⊥ , b ⊥ and r ⊥ comes from the fact that these transverse coordinates have the same physical meaning as those for the self-energy crossing the shockwave. After this change of variables, one can write the cross-section as with ω = z g /(z 2 (z 1 − z g )). The extra factor z 1 /(z 1 − z g ) in the second line comes from the Jacobian of the change of coordinates. The back-to-back limit of this expression can now be obtained as previously from the "correlation expansion" of the color correlator Ξ NLO,1 , keeping the leading terms in u ⊥ , u ⊥ and r ⊥ in the regime u ⊥ ∼ u ⊥ ∼ r ⊥ . We find for the leading term, This time however the term proportional to r i ⊥ does not vanish since the r ⊥ integral depends on the scalar product r ⊥ · u ⊥ and not only on r 2 ⊥ . Thanks to this expansion, the leading term can be written in a factorized form up to corrections of order q ⊥ /P ⊥ or Q s /P ⊥ : with the hard factor defined to be (3.55) and the integrals I 3 and I i 4 given by As for the case of SE 1 , we cannot make any approximations in these integrals since the correlation expansion of Ξ NLO,1 assumes r ⊥ ∼ u ⊥ . In particular, ∂ i P ⊥ I 3 and I i 4 are both "leading power" since they both fall like 1/P 3 ⊥ at large P ⊥ . This confirms the importance of the second term in the correlation expansion Eq. (3.53) of Ξ NLO,1 . On the other hand, any higher order corrections to Eq. (3.53) give a result vanishing at least like 1/P 4 ⊥ . Thus such contributions are suppressed in the back-to-back limit using the same rationale as in the previous subsection.
Remarkably, the two integrals I 3 and I i 4 can be computed analytically using tricks similar to those employed in the calculation of the hard factor in SE 1 . A detailed computation of these integrals is provided in Appendix A.4, and we only quote here the final results: (3.59) Plugging these expressions into Eq. (3.55), we obtain a representation of H NLO,4 as a single integral over z g .
Subtracting the kinematically constrained slow divergence. Similarly to SE 1 , this z g integral has both double and single logarithmic divergences at z g = 0, but is regularized by the counterterm H λ=L,ij NLO,4,slow (P ⊥ ) obtained by extracting the piece proportional to the WW gluon TMD in the subtraction term dσ λ=L V1,slow . The latter should include the kinematic constraint as discussed in the previous subsection: with the color correlator evaluated at the transverse coordinates As we have discussed, the Θ-function in Eq. (3.60) is the only modification as compared to the "naive" rapidity subtraction term giving leading-log BFKL/BK/JIMWLK evolution. Expanding around b ⊥ to extract the WW gluon TMD, we get the result for the subtraction term for slow gluons with the integrals I 3,s and I i 4,s given by (3.63) As outlined at the beginning of this subsection, it is crucial for the consistency of our kinematically improved evolution that the same choice of Q f as for SE 1 gives a welldefined z g integral for V 1 . The calculation of the two integrals I 3,s and I i 4,s is also detailed in Appendix A.4. They read, up to powers of z g corrections: Comparing these expressions with Eqs (3.58) and (3.59), one verifies that the kinematically constrained slow gluon subtraction term dσ λ=L V1,slow gives precisely the counterterm H λ=L,ij NLO,4,slow (P ⊥ ) that makes the z g integral in H λ=L,ij NLO,4 (P ⊥ ) convergent, provided we choose (3.67) The relation between z f and Q f can be seen to be identical to that derived from SE 1 in Eq. (3.48).
Final result. With the inclusion of the subtraction term for slow gluons H λ=L,ij NLO,4,slow (P ⊥ ) that we just computed, the remaining z g integral in Eq. (3.55) is convergent. The calculation of this integral is rather tedious. In the end, our final result for the hard factor associated with V 1 is, (with χ =Q/P ⊥ = Q/M qq ): The last two lines come from the slow gluon phase space between z f and z 1 (respectively z 2 for the quark-antiquark exchanged diagram V 1 ) which is not subtracted off since Eq. (3.62) only subtracts the phase space below the factorization scale z f . In particular, one notices that the C F term (left unintegrated in z g ) cancels against the corresponding term for diagram SE 1 (the second line in Eq.(3.49)) in such a way that the z f dependent piece of the sum of C F H NLO,3 − H NLO,4 is proportional to the color factor N c , in agreement with the overall color factor of the evolution equation of the WW gluon TMD. In the last line, one notices a double logarithmic dependence upon the factorization scale z f which will cancel in our final result with the ln 2 (Q 2 f /P 2 ⊥ ) term in Eq. (2.33) once Q f is expressed in terms of z f and the kinematics of the dijet process thanks to Eq. (3.67). The absence of ln 2 (z f ) dependence in the impact factor is another important feature of TMD factorization at NLO at small x Bj .

Soft contributions from real diagrams
In the previous section, we discussed the virtual contributions to the back-to-back NLO cross-section and demonstrated that they factorize into hard factors and the WW gluon TMD in the back-to-back limit P ⊥ q ⊥ , Q s . We explicitly computed each of the hard factors obtaining analytical results for them. Further, we showed that the leading log rapidity evolution of the WW TMD distribution must satisfy a kinematic constraint corresponding to lifetime ordering of slow gluons. Not least, we were able to obtain a unique transverse momentum factorization scale Q f that only depends on the external kinematic variables of the back-to-back cross-section and on the rapidity factorization scale z f . This section is dedicated to real contributions to the back-to-back NLO cross-section. We review the results from [1] for the contribution to the back-to-back cross-section from soft gluon radiation emitted from quark/antiquark after the interaction with the shockwave in the amplitude and the complex conjugate amplitude. Furthermore, we revisit the real emissions where the gluon interacts with the shockwave either in the amplitude or in the complex conjugate amplitude and we argue that they give power suppressed contributions in back-to-back kinematics.

Final state soft gluon radiation: Leading power Sudakov logarithms and O(α s ) corrections
Amongst all the real diagrams gathered in Fig. 2, only the subset with gluon emissions after the shockwave contribute to the back-to-back dijet cross-section at NLO and leading power in q ⊥ /P ⊥ . These contributions are labeled R 2 × R 2 and R 2 × R 2 , plus those obtained after quark-antiquark exchange. The diagram R 2 × R 2 contributes in two ways to the back-to-back dijet cross-section. Firstly, when the real gluon is emitted inside the quark jet cone, the phase space integration is divergent because of the QCD collinear singularity but this divergence cancels with the virtual cross-section. The finite leftover contribution is included in Eq. (2.33). It also contributes via soft gluons emitted outside the jet cone as [1] The R dependence of this result comes from the out-of-cone condition which obviously depends on the jet radius. Similarly, the diagram R 2 × R 2 contributes in the soft regime via the term Since this diagram has no collinear divergence, the in-cone phase space is power of R 2 suppressed and to leading order in R, the soft out-of-cone component does not depend on the jet radius. Both the expressions above display a dependence on the WW gluon TMD but the phase e −iξP ⊥ ·r bb and the logarithm couple the P ⊥ and r bb dependence. The ξ integral was calculated in [1]. It gives Sudakov double and single logarithms which are resummed in the Sudakov form factor S(P 2 ⊥ , r 2 bb ) in Eq. (2.33). It also contributes finite pieces to the coefficient function and these are included in Eq. (2.33). Soft gluon radiation has an important effect on azimuthal anisotropies because of the scalar product P ⊥ · r bb in the phase. Therefore the finite terms coming from soft gluons emitted after interacting with the shockwave, both in the amplitude and in the complex conjugate amplitude, depend on the harmonics of the Fourier decomposition of the dijet cross-section with respect to the azimuthal angle between P ⊥ and q ⊥ . This Fourier decomposition is given by dσ (n),λ (P ⊥ , q ⊥ , η 1 , η 2 ) cos(nφ) , (4.3) with φ the azimuthal angle between P ⊥ and q ⊥ . After absorbing the Sudakov logarithms into the Sudakov form factor, the finite pieces coming from Eq. (4.1) and Eq. (4.2) contributing to the coefficient function (and multiplying the linearly polarized WW distribution) take the form dσ (0),λ=L for the azimuthally averaged cross-section. These contributions, and their quark-antiquark exchange counterparts, are included in Eq. (2.33). For the cos(2φ) anisotropy, the corresponding results are dσ (2),λ=L with θ the angle between q ⊥ and r bb in polar coordinates. The finite terms for higher harmonics coming from soft gluon radiation have been computed in Appendix E of [1]. To sum up, real soft gluon radiation emitted after the shockwave contributes to the back-to-back dijet cross-section to leading power and their contribution is proportional to the WW gluon TMD. They also induce azimuthal anisotropies in such a way that the azimuthally averaged cross-section depends on the linearly polarized WW TMD and the cos(2φ) anistropy is sensitive to the unpolarized WW TMD. As shown in Section 5.1 of [1], the contribution from hard gluons with transverse momenta k g⊥ z g P ⊥ to diagrams R 2 × R 2 and R 2 × R 2 is power suppressed in the back-to-back limit.

Real emissions crossing the shockwave: Beyond leading power O(α s ) corrections
In order to prove factorization of the back-to-back dijet cross-section at NLO in terms of the WW gluon TMD, we must also demonstrate that the other real diagrams, in which the gluon crosses the shockwave either in the amplitude or in the complex conjugate amplitude are power suppressed in the back-to-back limit. The main difficulty with these diagrams comes from their dependence on CGC correlators which do not naturally collapse to the WW gluon TMD as the leading order one does in Eq. (2.11). This was also apparently the case for the virtual diagrams in which the gluon scatters off the shockwave; however as we discussed in Section 3, there exists a correlation expansion that enables one to isolate the leading power term which depends only on the WW gluon TMD.
This correlation expansion relies on the hierarchy of transverse scales u ⊥ ∼ r ⊥ b ⊥ , where typically u ⊥ ∼ 1/P ⊥ , r ⊥ ∼ 1/|k g⊥ | and b ⊥ ∼ 1/q ⊥ . In particular, for virtual graphs with gluons crossing the shockwave, this expansion gives a mathematically well-defined impact factor thanks to the kinematic constraint. Effectively, the kinematic constraint enables one to clearly distinguish gluons with hard momenta, of order of P ⊥ , that contribute to the impact factor, from gluons with momenta of order Q s , q ⊥ , that control the WW gluon TMD and its evolution at small-x. Hence, the contribution to the NLO back-toback impact factor of virtual diagrams with gluons crossing the shockwave is dominated by hard gluons with transverse momentum of order of P ⊥ (such gluons are not kinematically forbidden in back-to-back kinematics since the gluon is virtual). For real emissions, one can similarly expand the color correlator. Imposing back-to-back kinematics requires |k g⊥ | to be smaller than q ⊥ , or in terms of conjugate transverse coordinates, r ⊥ larger than b ⊥ . One may then question the validity of the correlation expansion. However, the kinematic constraint, which is naturally enforced in real graphs with a gluon crossing the shockwave through the argument of the Bessel function, enforces that |k g⊥ | √ z g P ⊥ q ⊥ , i.e. r ⊥ b ⊥ . Therefore, the correlation expansion remains well defined but the lack of overlapping phase space between real gluons satisfying both the kinematic and back-to-back constraints leads to power suppression of these diagrams 7 . In Appendix B.2, we provide a more mathematically grounded demonstration showing that diagrams R 1 × R 2 , R 1 × R 2 , R 1 × R 1 and R 1 × R 1 are indeed power suppressed in the back-to-back limit.

Weizsäcker-Williams gluon TMD factorization at NLO from the CGC
In this section, we will put together all of the results of the previous sections for the NLO back-to-back dijet cross-section that we have shown can be expressed as remarkably simple factorized expressions for the convolution of the NLO impact factor and the WW TMD gluon distribution, with the former expressed as a product of a universal soft factor and a coefficient function. The WW distribution obeys a kinematically constrained JIMWLK rapidity evolution equation.
In Section 5.1, we will discuss the structure of these results, and outline subtle features that emerge when one considers the rapidity evolution of the WW TMD distribution. In addition, we show that one needs to further reorganize our results using an appropriate choice of a hard factorization scale of the Sudakov form factor to obtain physically sensible cross-sections. Our final NLO results for all moments of the back-to-back dijet cross-section are compactly summarized in Section 5.2.

Structure and reorganization of TMD factorized results
We first combine the analytic results of the previous sections for the functions H NLO,1 and H NLO,2 with those that as we have seen greatly simplify the "other" contribution dσ (0),λ=L other in Eq. (2.33) into two new hard factors H NLO,3 and H NLO,4 . For the purpose for a discussion of the choice of the hard scale, we shall first write the fixed order O(α s ) result without assuming Sudakov exponentiation; we will expand the Sudakov form factor in Eq. (2.33) to first order in α s . For jets defined with a generalized k t algorithm, we find that our previous result in Eq. (2.33) becomes , given by for a longitudinally polarized virtual photon. These two functions also depend on z 1 , z 2 and R but we choose here to leave manifest only the dependence upon the dimensionless ratio Q/M qq and the rapidity factorization scale z f . Because of theQ 2 factor included in H 0,λ=L LO (P ⊥ ), the photo-production limit Q 2 → 0 of the cross-section is zero although the functions f λ=L 1 and f λ=L 2 diverge logarithmically as χ → 0: Albeit finite, the χ → 0 limit or equivalently the M 2 qq Q 2 limit, leads to an unphysical negative NLO cross-section when α s | ln(χ)| ∼ 1. This unphysical behaviour, resulting from imposing an additional hierarchy of scales M 2 qq Q 2 on top of the W 2 Q 2 and P ⊥ q ⊥ conditions, points towards the need for a more sytematic study of higher order corrections or some additional resummation of ln(χ) logarithms. This will be also discussed below in the opposite regime Q 2 M 2 qq . For the record, and for later convenience, we can also define the NLO hard factor associated with virtual gluon emission and unresolved real emissions (those inside the jet cone): Rapidity factorization scale z f dependence. The NLO coefficient function has a remaining rapidity factorization scale z f dependence. Remarkably, this dependence enters through the z f term in the finite piece only, without double logarithms of z f nor Q f . Indeed, the transverse momentum scale Q f introduced in the kinematically constrained rapidity evolution of the WW gluon TMD has been simplified thanks to the relation As noted previously, this constraint follows from our significantly more detailed analysis here of the virtual diagrams in the back-to-back limit, and should be employed when evolving the TMDsĜ ij Y f from some initial scale Y 0 = ln(z 0 ) up to the factorization scale Y f = ln(z f ) with the evolution equation given by Eq. (2.37). Eq. (5.7) also enables one to simplify the single logarithmic coefficient s f in the Sudakov form factor, defined in Eq. (2.36). Plugging the identity Eq. (5.7) into s f , we find and note thus that s f does not depend anymore on the rapidity factorization scale z f . Our study enables us to unambiguously fix the value of the Sudakov logarithms to single logarithmic accuracy. Eq. (5.7) is also consistent in the photo-production limit since z f /Q 2 f then goes to ec 0 /M 2 qq . With regard to the z f dependence of the result, one might be tempted to derive the Y f = ln(z f ) evolution equation of the WW gluon TMD by requiring that the cross-section should be independent of Y f up to O(α 2 s ) sub-leading corrections. Then from Eq. (5.1), one would obtain the linear RG equation for the WW gluon TMD: This is however not the particular evolution equation Eq. (2.37) mentioned in the first section that was derived employing the JIMWLK equation in [95] (see also [98]). Since we noted there that this equation does not generate a closed form expression, clearly something is missing here. The resolution lies on the observation that the kinematic constraint interpolates between full non-linear evolution in Eq. (2.37) when z g 1 to the linear evolution in Eq. (5.9) when z g ∼ O(1). As we will see below, this is because the correlation limit expansion is only valid in the NLO impact factor where z g ∼ O(1), whereas it breaks down in the small-x enhanced terms where z g 1. In deriving our impact factor, we implicitly assume z f ∼ O(1); thus all the non-linearities are absorbed into the small-x evolution ensuring the impact factor factorizes from the WW distribution.
To illustrate this, we can formally write the fixed order NLO cross-section as a sum of the LO cross-section, the NLO impact factor and the O(α s ln(x 0 /x f )) slow contributions: c.,slow (z g ) , where I k.c.,slow incorporates the slow gluon divergence including the kinematic constraint. The second term in the expression on the r.h.s corresponds to the impact factor. Note that I(z g ) implicitly contains step functions that bound the z g integral from above. The third term on the r.h.s is the large ln(z f /z 0 ) ∼ ln(x 0 /x f ) rapidity logarithm. At this step, the cross-section is indeed obviously independent of z f . Now we can schematically write the correlation expansion as I(z g ) = I (0) (z g ) + I (1) (z g ) + ... + I (n) (z g ) + ... , (5.11) k.c,slow (z g ) + I (1) k.c,slow (z g ) + ... + I (n) k.c,slow (z g ) + ... , (5.12) where I (0) represents the leading power contribution proportional to the WW gluon TMD and the terms I (n) are suppressed by (q ⊥ /P ⊥ ) n . By construction, the kinematically constrained slow gluon counterterms I (n) k.c,slow (z g ) cancel the z g → 0 divergence order by order in the correlation expansion. Hence in the impact factor term, the leading power term in the back-to-back limit is simply This is the leading contribution in the correlation expansion computed analytically in this paper from all virtual diagrams and from real soft emissions emitted after the shockwave interaction with the quark-antiquark dipole. However the correlation expansion breaks down for the third term in the r.h.s of Eq. (5.10) when z 0 is small and all I (i) k.c,slow terms contribute equally. This is because the correlation expansion relies on r ⊥ ∼ u ⊥ b ⊥ , b ⊥ , with each additional term in Eq. (5.12) having an extra power of r ⊥ or u ⊥ . Therefore I (n) k.c,slow (z g ) behaves parametrically like (for n ≥ 1) . Since Q f ∼ P ⊥ , one may be tempted to say that the larger n is, the more it is suppressed by powers of P ⊥ . Yet the kinematic constraint, represented by the step function which couples the r ⊥ and z g dependence, makes each term more divergent in the z g → 0 limit. One is thus forced to keep the full I k.c,slow integrand in the phase space between z 0 and z f , leading to The first two terms on the r.h.s correspond to the result Eq. (5.1), while the last one gives the first step of the rapidity evolution of the WW gluon TMD. It is now clear that the full cross-section, as written, is not strictly z f independent. Hence the evolution equation for the TMD in the LO cross-section should be inferred from the form of the last term in Eq. (5.17). Since this term is precisely the first step of the kinematically constrained JIMWLK evolution equation for the WW distribution (Eq. (2.37)), as demonstrated in [1,37], the all order resummation of small-x logarithms should be performed using this evolution equation. If one insists on writing the impact factor such that one gets Eq. (2.37) from the z f independence of the NLO cross-section up to O(α 2 s ) corrections (as one must), we need to add the to the result in Eq. (5.1), the term Although this term is superficially power suppressed in the back-to-back limit for z f ∼ 1 (as it should be), it allows us to recover the non-linear evolution of the WW gluon TMD Eq. (2.37) by taking the Y f = ln(z f ) derivative of the NLO cross-section. When z 0 is not small (x Bj ∼ 1), it makes sense to truncate the correlation expansion of I k.c.,slow to first order. In this case, the Y f dependence of the WW gluon TMD should be given by the linear equation Eq. (5.9). However this evolution equation (given in Eq. (5.9)) is not meaningful since when x Bj ∼ 1 there is no need for small x resummation in the first place.
Jet definition. The NLO coefficient function displays a ln(R) dependence which is the remnant of the infrared and collinear safe jet definition of the final state. The formula in Eq. (5.1) is valid for jets defined using the generalized k t algorithm. (These include the anti-k t algorithm [99] and the Cambridge/Aachen algorithm [100,101].) For such jet-cone algorithms, one should subtract from the NLO hard factor the quantity [102] We leave for future phenomenological study the computation of the NLO coefficient function and soft factor for dihadrons instead of dijets. In that case, the NLO coefficient function would not depend on R but rather on the factorization scale µ F corresponding to the evolution of the dihadron fragmentation function [43].
The regime Q 2 M 2 qq . Besides Sudakov logarithms which are resummed via the Sudakov form factor, the impact factor may potentially contain large logarithms that spoil the perturbative convergence of the cross-section or lead to unphysical results such as negative cross-sections. One example is the ln(R) term which becomes large for jets defined with small jet radii. The resummation of such logs can be performed using "micro-jet" evolution [103] or jet functions in the Soft Collinear Effective Field Theory (SCET) approach [102].
Another potentially problematic logarithm involves the terms containing ln(χ) (with χ = Q/M qq ), when χ becomes large. Indeed, χ can in principle be arbitrarily large for a given P 2 ⊥ , while maintaining the back-to-back constraint 8 . In this χ 1 limit, the NLO hard factor behaves as (5.20) When χ 1, the NLO cross-section becomes negative due to the dominant −4 ln 2 (χ) term in the impact factor. Such negative double logarithms of the ratio between the hard scale in the Sudakov logarithms (here P ⊥ ) and the photon virtuality Q also appear in lepton-jet correlation in DIS [104] or in back-to-back heavy quark pair production in DIS [105] within the TMD framework.
In our calculation, the origin of this large − ln 2 (χ) contribution to the NLO crosssection can be traced to our choice of the hard scale µ 2 h in the Sudakov logarithms in Eq. (5.1) to be P 2 ⊥ . Indeed, the double logarithm is a leftover of the kinematically improved small-x evolution given by the last term of the third line in Eq. (2.33), namely Since Q 2 f ∼ Q 2 when χ 1, this is the negative double logarithm we were looking for. The presence of P 2 ⊥ in the denominator of the argument of this double logarithm comes from the choice of P 2 ⊥ as the hard scale in the Sudakov double logarithm. 8 This would correspond to the kinematic domain P 2 Having identified the origin of this term, we note that we have the freedom to make a different choice for the hard scale in the Sudakov evolution to ensure that the cross-section is positive definite and the coefficient function remains of order O(α s ) for this scale choice. We will therefore change the hard scale in the Sudakov resummation to instead of P 2 ⊥ . In doing so, we ensure that the contribution to the large logarithm from the phase spaceQ 2 ∼ µ 2 h P 2 ⊥ is resummed in the Sudakov factor. Introducing this choice of hard scale µ h inside the Sudakov double and single logarithms of the first line in Eq. (5.1), we can rewrite the NLO O(α s ) Sudakov logarithms We can now put the terms in the second line above in the impact factor while the terms in the first line can be exponentiated into the Sudakov form factor as with s 0 = ln 2(1 + cosh(∆η 12 )) The modified hard factor can similarly be expressed in terms of the functionsf 1 and f 2 : With this modified Sudakov form factor, the behaviour of the NLO hard factor at large χ 1 becomes where the coefficient of the double logarithm ln 2 (χ) is now positive. Hence by changing the hard scale of the Sudakov form factor, we were able to obtain a cross-section that is positive definite even in the regime χ 1 (Q 2 M 2 qq ). For an estimate of the effect of this change of scale, we refer the reader to Fig. 8 in Section 6.
Although positive definite in the regime χ 1, the NLO coefficient function still has double and single logarithms of χ which may spoil the perturbative convergence of the cross-section when enforcing the additional constraint Q 2 M 2 qq . A two loop calculation of the NNLO coefficient function would help to clarify this issue, or even better, stabilize the cross-section as observed in other processes [106,107].
Anomalous dimension of the WW gluon TMD. Besides the NLO diagrams computed in sections 3 and 4 which include all O(α s ) quantum corrections obtained using the small fluctuation propagator (CGC effective vertices) in the classical background field, there is an additional contribution in the CGC EFT coming from the one-loop correction to the classical background field A +,a cl itself [92,[108][109][110]. The diagram associated with this contribution is displayed in Fig. 4. As shown in [92], this diagram has a UV divergence, which is the same as the standard UV divergence of the gluon vacuum polarization diagram at one-loop, namelyĜ  This UV divergence is simply removed by replacing the bare coupling constant in the definition Eq. (2.23) of the WW gluon TMD by the renormalized one [92]. The TMD then acquires an additional scale dependence µ given bŷ One may then ask what the natural scale µ should be. The answer to this question depends on the hard scale which is probed by the process to be considered. At the level of the WW gluon TMD alone, µ should be set by the saturation scale or q ⊥ . Since the WW gluon TMD only depends on r bb (assuming a homogeneous target in the transverse plane), a reasonable choice would be µ 2 = c 2 0 /r 2 bb , where the c 0 constant is rather arbitrary here. Such ambiguities can be absorbed into the initial condition for the WW gluon TMD.
However, for the back-to-back hard dijet process, the natural hard scale is rather µ h ∼ P ⊥ , and therefore, α s in the NLO corrections should run at this specific value. The phase space between the "initial" scale c 2 0 /r 2 bb of the TMD and the hard scale of the process leads to an additional Sudakov single logarithm, as we now demonstrate. The one-loop running coupling satisfies Plugging this expression for the WW gluon TMD renormalized at the hard scale µ h into the LO cross-section Eq. (2.25) gives a NLO correction to the impact factor, which displays a single Sudakov logarithm with coefficient β 0 . This Sudakov single log comes physically from the large phase space for the logarithmic running of the coupling between P ⊥ and q ⊥ ∼ 1/r bb . From the one-loop QCD β-function and Eq. (5.31), one can readily deduce the µ evolution equation that enables one to resum this single logarithmic correction: Integrating this equation between the natural initial evolution scale c 2 0 /r 2 bb of the TMD and the hard scale µ 2 h results in the following exponential which can be included into the Sudakov form factor as it provides a single logarithmic correction of the same order as the s 0 or s f term in Eq. (5.25). The somewhat arbitrary choices of factorization scale µ h and initial scale c 2 0 /r 2 bb is equivalent to a choice of scheme to remove the UV divergence of the diagram in Fig. 4. In practice, one can gauge the scheme dependence by varying the renormalization scale µ h by a factor of two around our central choice µ 2 h = P 2 ⊥ +Q 2 in our final expressions given in section 5.2. Including in Eq. (5.25) the anomalous dimension of the WW gluon TMD Eq. (5.35) leads to the following Sudakov form factor: The presence of the β 0 term is in agreement with the calculation done in [111] within the CSS formalism. However our result differs from the Sudakov factor obtained in conventional collinear factorization [78,80] by the 1 term in the single log coefficients f . This difference (which arises in our case from the kinematically constrained small-x evolution for the WW gluon TMD) deserves further investigation.

Final results for moments of the NLO back-to-back dijet cross-section
Following the extended discussion in Section 5.1, and the subsequent reorganization of our results, we can now finally write down the principal results of this paper for the moments of the NLO back-to-back longitudinally polarized dijet cross-section. For the zeroth moment of the azimuthal anisotropy, we obtain The functionf λ=L 1 is given by Eqs. (5.2) and (5.27) while andf λ=L 2 is given by (5.3) and (5.28). The Sudakov form factor is given by Eq. (5.36). In writing this expression, we exponentiated our fixed order results for the Sudakov double and single logarithms to all orders presuming that CSS evolution can be extended to small x. A formal proof of this exponentiation to all orders does not exist at small x and is an important topic for future research.
In a similar fashion, by collecting the various terms in Eq. (2.34) computed in Section 3, we obtain the analytic expression (modulo the Fourier transform of the WW gluon TMDs) for the cos(2φ) anisotropy: For completeness, we recall here the results obtained in [1] for higher even harmonic coefficients. They arise solely from real soft gluon radiation as virtual corrections do not generate 9 cos(nφ) anisotropies for n ≥ 4. They read where h(p) = p k=1 1 k is the p th harmonic number. Our results significantly extend those obtained previously in [76][77][78][79][80]. Results for the back-to-back dijet cross-section for transversely polarized virtual photons will be presented in a separate work.
6 Numerical analysis of the results in Section 5.2 In this section, we will provide a brief numerical study of the results summarized in Section 5.2. Our aim here is to illustrate qualitatively the relative magnitude of Sudakov suppression and the NLO coefficient function to the impact factor. In our numerical study we will employ the McLerran-Venugopalan (MV) model to compute the Weizsäcker-Williams (WW) gluon distribution. We will present a more detailed examination of the interplay of Sudakov suppression and the energy evolution of the WW gluon distribution in a subsequent publication.
For simplicity, we shall neglect the impact parameter dependence of the nuclear target; the differential cross-section is therefore proportional to the effective area S ⊥ of the target. We will consider here the differential yield which we can decompose into Fourier modes with respect to φ (the angle between P ⊥ and q ⊥ ) as We will present results for the azimuthally averaged differential yield dN (0),λ=L , the elliptic anisotropy v λ=L 2 , and the quadrangular anisotropy v λ=L 4 . Before presenting our numerical results we will first discuss the nonperturbative modeling of the WW gluon TMD at small-x and nonperturbative contributions to the Sudakov form factor.

The WW gluon TMD in the MV model
While the rapidity evolution of the WW TMD satisfies the kinematically constrained JIMWLK renormalization group equation, the initial conditions for the evolution are intrinsically nonperturbative. To compute the WW gluon TMD, we will rely on a Gaussian approximation which allows one to express any correlator of Wilson lines (or their derivatives) in terms of a nontrivial funciton of the two point dipole correlator of Wilson lines S Y (r bb ) [112][113][114]. Within this approximation, the WW gluon TMDĜ ij Y (r bb ) can be expressed as [11,19,25] with Γ(r bb ) = − ln (S(r bb )) .

(6.4)
For brevity, we dropped the rapidity dependence Y in the dipole and the WW correlators as we only consider the MV model in our numerical analysis. In our numerical study, we will use the MV model parametrization where we choose Q 2 s0 = 1.0 GeV 2 and a nonperturbative scale m = 0.241 GeV of order Λ QCD (see also [115,116]) as a proxy for a large nucleus.

Modeling the Sudakov form factor
The Sudakov form factor is computed using Eq. (5.36) with the two loop running coupling. Our default choice for the hard scale is µ 2 h = P 2 ⊥ +Q 2 . The analytic expressions are provided in appendix C. When computing the Sudakov factor with running coupling, one has to make sure that the QCD Landau pole does not belong to the µ integration range. In other words, we will ensure that . (6.6) In the TMD literature, there is a widely used prescription to tame the nonperturbative contributions from large distances [70]. The lower bound of the µ 2 integration in Eq. (5.36) is changed from c 2 0 /r 2 bb to c 2 0 /r * 2 bb where In order to compensate for the missing contribution at large r 2 bb , one includes a nonperturbative Sudakov contribution following [117,118]  with the overall factor of 2 in the exponent due to the two quark jets in the final state. We use the same values as in [117,118] (see also [119,120] in the context of small-x physics) for g 1 = 0.212 GeV 2 , g 2 = 0.84, Q 2 0 = 2.4 GeV 2 , and r max = 1.5 GeV −1 , which are values constrained from fits to experimental data.

Numerical results
The NLO impact factor has two contributions: (i) the soft Sudakov factor containing double and single Sudakov logs, and (ii) the O(α s ) coefficient function which includes the NLO hard factor as well as finite pieces from soft gluon emissions outside the jet cone.
We begin this section by first studying the impact of the Sudakov form factor on the azimuthally averaged cross-section and the elliptic anisotropy as a function of q ⊥ . Previous phenomenological studies for dijet/dihadron production in DIS at small x [73][74][75] only considered the Sudakov double logs for fixed strong coupling and ignored the nonperturbative contribution to the Sudakov form factor. In Fig. 5, we will illustrate the effect of Sudakov suppression beyond the double log approximation. For this purpose, we shall examine the LO impact factor + contributions to the NLO impact factor containing Sudakov logarithms (and excluding the NLO coefficient function): and dσ (2),λ=L = α em α s e 2 f δ (2) We plot Eqs. (6.10) and (6.11) with four different contributions to S(µ 2 h , r 2 bb ): (i) double log with fixed coupling, (ii) double and single logs with fixed coupling, (iii) double and single logs with running coupling, and (iv) double and single logs with running coupling and with non-perturbative Sudakov factor. For the fixed coupling case, we evaluate the coupling at the hard scale µ h . For reference, we also include the LO impact factor which is obtained by setting the Sudakov form factor to 1 in Eqs. (6.10) and (6.11). Firstly, as also noticed in [73][74][75], the Sudakov form factor including double logarithms alone leads to a significant suppression of the yield at small q ⊥ . For out choice of kinematics and jet radii, adding the contribution of the Sudakov single log results in enhancement of the cross-section and its elliptic anisotropy as its coefficient has the opposite sign to that of the double log 10 . In addition, we observe a significant suppression of the cross-section when we switch from fixed coupling α s (µ 2 h ) to the running coupling (at two loops) α s (µ 2 ). Since α s (µ 2 h ) < α s (µ 2 ) when µ 2 h > µ 2 , soft gluon emission is enhanced in the running coupling case and results in a stronger Sudakov suppression. The effect of the running of the coupling can be understood in more detail by expanding Eqs. (C.8): where L = ln . The fixed coupling case accounts only for the first terms in Eqs. (6.12) and (6.13). On the other hand, the running of the coupling at two loops contains all leading logs (LL) α n s L n+1 , next-to-leading-logs (NLL) α n s L n , as well as a subset of the N 2 LL terms. These contributions are responsible for the stronger Sudakov suppression observed in Fig. 5. Lastly, adding the non-perturbative contribution to the Sudakov factor further suppresses the cross-section. Since the effects plotted are large, to understand the interplay between saturation and Sudakov suppression it is therefore crucial that both double and single logs as well as the running of the coupling and the non-perturbative contribution to the Sudakov form factor are taken into account in phenomenological studies. In Fig. 6, we turn our attention to the O(α s ) NLO coefficient function that multiplies the soft factor. We compute the azimuthally averaged yield as a function of q ⊥ (in the left panel) and as a function of P ⊥ (in the right panel). We display the results employing the LO impact factor, the LO + NLO contributions to the impact factor that contain Sudakov logarithms (see Eq. (6.10)), and the LO + full NLO impact factor which contains both Sudakov and the finite O(α s ) coefficient function in Eq. (5.37). The size of the NLO coefficient function are strongly dependent on P ⊥ or more precisely, on the ratio χ =Q/P ⊥ as expected from the factorization formula in Eq. (5.37). At small values of χ the NLO finite pieces enhance the differential yield while they decrease it at large χ. As noted in [1], the finite NLO terms for the azimuthally averaged cross-section have a contribution from linearly polarized gluons due to their interplay with soft gluon radiation. This contribution is rather small in the saturation region q ⊥ Q s since the relative contribution from linearly polarized gluons in the WW distribution is suppressed. At the end of day, the NLO coefficient function is nearly independent of q ⊥ and leads to a constant shift of the crosssection as a function of q ⊥ as shown on Fig. 6-left when going from the blue to the red curve.  Next we study the elliptic and quadrangular anisotropies as a function of the momentum imbalance q ⊥ in Fig. 7. At LO, as seen from Eq. (2.30), the elliptic anisotropy is proportional to the ratio between linearly polarized to unpolarized gluons which grows with q ⊥ . The Sudakov factor significantly decreases the elliptic anisotropy. This effect is nontrivial since v 2 is defined as a ratio and one expects therefore a partial cancellation of the Sudakov suppression. However the unpolarized and linearly polarized WW gluon TMDs are not identical and thus the Sudakov factor affects them differently. On the other hand, the effect of the NLO coefficient function is to enhance the elliptic anisotropy, as it now acquires a contribution from unpolarized gluons. This can be seen from the last two lines of Eq. (5.38) and physically corresponds to the effect of soft gluon emissions that are near to but yet outside the jet cone. These emissions tend to preferentially align the imbalance q ⊥ parallel (or anti-parallel) to the transverse jet direction P ⊥ [80] (see also [121] for the case of lepton-jet correlations). The NLO finite pieces also produce a small quadrangular anisotropy (5.40) which is completely absent in the leading order and in the LO + "Sudakov only" result. We close this section by studying the Q dependence of the azimuthally averaged differential yield for two different choices of the hard factorization scale µ h in the Sudakov factor. As discussed in Section 5.1, different choices of µ h also lead to different NLO coefficient function. A common choice of the factorization scale is µ 2 h = P 2 ⊥ [73,77,122]. However, as we discussed in Section 5.1, we find that the associated NLO finite piece develops large negative logs at large Q, resulting in an unphysical negative yield. A more natural choice is the combination µ 2 h = P 2 ⊥ +Q 2 which renders the yield positive. As expected, these scheme choices do not affect the Q P ⊥ region.

Summary and outlook
In this paper, we performed the first complete analysis of the back-to-back limit of the NLO dijet impact factor at small Bjorken x Bj in the CGC effective field theory. Our study significantly extends the results obtained in our previous paper [1] where the emphasis was placed on Sudakov double and single logarithms appearing in the NLO impact factor and the interplay of these contributions with rapidity evolution. It was believed in [1] that NLO hard factor contributions not containing Sudakov logarithms were only factorizable in the strict dilute TMD limit of P ⊥ q ⊥ Q s . Our principal result here is that the back-to-back limit of the NLO impact factor can indeed be written in a fully factorized form to all orders in Q s /q ⊥ and up to power corrections in q ⊥ /P ⊥ (or equivalently, Q s /P ⊥ ).
This factorized expression involves three important ingredients. (i) The WW gluon TMD which, in light cone gauge and light cone quantization, represents the number density of gluons per unit momentum in the target at small-x. This WW TMD obeys a kinematically constrained JIMWLK renormalization group equation in rapidity.
(ii) A Sudakov soft factor that resums to all orders double and single logarithms of P ⊥ /q ⊥ appearing in the impact factor at NLO. (iii) A NLO coefficient function that gathers all the O(α s ) finite terms which are not enhanced by rapidity logarithms or logarithms in P ⊥ /q ⊥ .
Following [1], we confirmed from a detailed analysis of all virtual diagrams that the WW gluon TMD must satisfy a kinematically constrained JIMWLK evolution equation in rapidity. We were able to go further and unambiguously fix the transverse scale in the kinematic constraint that imposes lifetime ordering of successive gluon emissions. We leave for forthcoming work a detailed numerical study of Eq. (2.37) giving the kinematically constrained evolution of the WW gluon TMD.
Fixing the transverse scale of the kinematic constraint also enables us to eliminate ambiguities in the single Sudakov logarithms. In the end, both Sudakov double and single logarithms turn out to be independent of the small-x resummation scheme. This extends to single logarithmic accuracy the results in [77] in which the authors claim that Sudakov and small-x resummation can be performed independently and simultaneously in the Regge limit.
We computed here for the first time remarkably compact analytical expressions for the NLO coefficient function for longitudinally polarized virtual photons. Such compact results will be very beneficial for quantitative predictions of our NLO results for future data from the Electron-Ion Collider. Our NLO results however display a residual rapidity factorization scheme dependence that can be quantified further beyond NLO.
We performed a preliminary numerical study of our analytic results with the aim of providing the reader with a feel for the order of magnitude of individual contributions in our factorized formula at NLO coming from the soft and hard contributions to the impact factor. We demonstrated the significant impact of the Sudakov double logarithmic soft factor contributions relative to the LO impact factor and further showed the relative impact of the contributions from the Sudakov single logarithms to the soft factor. We observed that the running of the coupling in the soft factor has a sizable impact on the differential cross-sections as a function of the dijet imbalance q ⊥ .
The effect of the NLO coefficient function can be summarized by distinguishing between the NLO hard factors coming from hard virtual and in-cone real gluons and from the rest of the contributions to the coefficient function that are associated with real soft outof-cone gluon emissions. On the one hand, the NLO hard factors do not modify the shape of the cross-section as a function of q ⊥ (in sharp contrast to the Sudakov form factor), and contribute an overall shift which can be comparable in magnitude to the soft factors. However, we observed that these contributions can qualitatively modify the shape of distributions as a function of P ⊥ since the NLO hard factor depends on the ratio P ⊥ /Q. On the other hand, finite pieces contributing to the NLO coefficient function that come from real soft gluons emissions outside the jet cone significantly alter the shape (as a function of q ⊥ ) of the v 2 and v 4 azimuthal anisotropies. A fully quantitative analysis of the contribution of these NLO finite pieces requires small-x evolution in order to mitigate the residual dependence of the hard factors on the rapidity factorization scale.
Our study can be developed further in several directions. Firstly, we can compute results for the hard factors in the scattering of transversely polarized virtual photons off the target to complement the results presented here for longitundinally polarized virtual photons. Secondly, it is crucial to compute the kinematically constrained small-x evolution of the WW gluon TMD extending prior work on the dipole distribution [34,123]. It is also important to understand the contribution of residual logarithms of the ratio of the virtuality of the photon and the dijet invariant mass, which can in principle be large for a wide separation of the two scales. This comment also applies for the remaining ln(R) dependence in the impact factor; for such contributions, we note that resummation techniques are already available [102,103]. Another open question with regard to back-toback dijets is the magnitude of the higher twist q ⊥ /P ⊥ contributions that require one to go beyond the correlation limit and are not considered here 11 . In principle, such contributions can be systematically computed following the method outlined in Appendix B, resulting in the appearance of distributions beyond the WW TMD.
Not least, a related and interesting follow-up calculation relevant for EIC phenomenology would be to replace our dijet computation with one for dihadrons using fragmentation functions. Forthcoming papers will address the aforementioned important questions towards the goal of making quantitative predictions for the EIC with the aim of isolating the dynamics arising from the energy evolution of saturated gluons from the physics governing the suppression of radiation at the edges of phase space.
The building block of H NLO,1 is the integral where, in the second line, we have performed the angular integration and the change of variable u = P ⊥ u ⊥ . Using the standard "replica" trick we can reexpress the u integral in terms of the ordinary hypergeometric function 2 F 1 (a, b, c; z): To obtain the second line, we used results from section 13.45 in [127] for Weber-Schafheitlin type integrals. Using the identity [128] z z − 1 , (A.4) and 2 F 1 (a, 0, 1; z) = 1, the derivative with respect to b gives Then using 2 F 1 (a, b, b; z) = (1 − z) −a , the derivative with respect to a can be easily computed: The integral I 0 then reads The first NLO hard factor can be related to the integral I 0 since . (A.8) After taking the derivative of I 0 with respect to P ⊥ , one ends up with the same result as Eq. (3.20) obtained using a different method.

A.2 Computation of H NLO,2
This appendix deals with the calculation of the integral I i 1 defined in Eq. (3.8) by The key manipulation consists of changing the variable l ⊥ to l ⊥ − K ⊥ , where in the definition of J in Eq. (3.5) so that the phases within cancel except for a remaining e il ⊥ ·u ⊥ factor. One can then switch the order of integration to perform the u ⊥ integral first: (A.11) The remaining integral over l ⊥ can be performed in polar coordinates: .
(A.12) The integral over θ is interesting. First one notices that for = K ⊥ , it gives whereas for > K ⊥ , one gets 2π 0 dθ 2π cos 2 (θ) + K ⊥ cos(θ) Hence the θ integration leads to a discontinuity at = K ⊥ and the integrand vanishes over a set of zero measure. The condition = K ⊥ corresponds to the condition |k g⊥ − P ⊥ | = (1 − z g /z 1 )|P ⊥ | in terms of the virtual gluon transverse momentum k g⊥ [37]. In particular, the kinematic configuration z 1 k g⊥ − z g P ⊥ = 0, namely, a virtual gluon collinear to the quark jet, does not contribute to the integral for this specific diagram.
Based on these results, one must separate the integral into two pieces I i 1 = I i 1,< + I i 1,> , (A. 16) with (A.17) (A. 18) In I i 1,< we have removed the −i prescription since |K ⊥ | < ∆ V3 and therefore the support of the integral does not encounter the pole. This integral reads .

(A.19)
For I i 1,> , we use the formula to properly account for the −i prescription in the contour integration, and we find in the end Adding these two terms together, one finds Eq. (3.10).

A.3 Computation of H NLO,3
In this section, we will provide details on the calculation of the integral I 2 defined by Eq. (3.35): We first consider the second (UV subtraction) term, whose dependence on u ⊥ through the exponential prefactor complicates the calculation. It is simpler to compute the integral with a UV subtraction term proportional to with ξ > 0 a fixed number independent of u ⊥ with the idea being that one adds and subtracts this term inside the r ⊥ integral. (See also [31,37].) Since [129]  we can write the hard factor as The integral over u ⊥ in the second line is done using the identity in Eq. (D.3). To compute the integral over u ⊥ and r ⊥ of the first line in Eq. (A.25), the trick is to use the identity [1] , (A. 26) which introduces two supplemental transverse integrals with the benefit of putting the transverse coordinate dependence on u ⊥ and r ⊥ inside pure phases. The integration over u ⊥ gives a Dirac delta function that enforces l 1⊥ = P ⊥ so that one easily performs the u ⊥ and l 1⊥ integrals, leading to with ∆ 2 = ω(P 2 ⊥ +Q 2 ). In the integral over r ⊥ , the ξ dependent UV subtraction term still regulates the UV r ⊥ → 0 behaviour of the diagram, since at small r ⊥ , K 1 (∆r ⊥ ) ∼ (∆r ⊥ ) −1 . In the IR regime, when r ⊥ is large, the modified Bessel function kills the power divergence beyond the scale ∼ 1/∆. Since ∆ 2 ∝ ω ∝ z g , one sees that in the slow gluon limit, the r ⊥ integral is more and more sensitive to this IR divergent domain. This explains the double logarithmic divergence of the resulting z g integral.
The last step of the transverse integration consists of computing the r ⊥ integral. Using the result Eq. (D.5), we end up with As expected, the ξ variable cancels between the two terms since the ξ dependent term was added and subtracted in the first place. The derivative w.r.t P ⊥ is straightforward and leads to Eq. (3.37).

A.4 Computation of H NLO,4
In this appendix, we shall provide details of the computation of the two integrals I 3 and I i There is a slight abuse of notation in Eq. (B.18). Indeed although the hard factor H λ=L,ijk R 1,1 ×R 2 depends on ξ ≡ z g /z 1 through the ω variable, it is sufficient to consider the leading term in the z g → 0 limit since the hard factor factors out of the ξ integral and does not depend on ξ.
The form of the kernel in Eq. (B.16) brings in the additional r k bb /r 2 bb dependence inside the Fourier transform of the WW gluon TMD in coordinate space in Eq. (B.18). The two terms in the hard factor H λ=L,ijk R 1 ×R 2 (inside the curly bracket) come respectively from diagram R 1 × R 2 (real emission from the quark before the shockwave with absorption in the c.c. by the quark after the shockwave) and diagram R 1 × R 2 (real emission from the anti-quark before the shockwave with absorption in the c.c. by the quark after the shockwave). It is very important to consider these two diagrams together as each individual diagram has a power divergence ∝ 1/ω in the limit ω ∝ z g → 0 which is potentially troublesome. This is because such a power divergence combined with the phase e −iξP ⊥ ·r bb would lead to a P ⊥ /q ⊥ enhancement of the cross-section. Fortunately the power divergence cancels amongst the two diagrams with the net result up to powers of z g corrections which are sub-leading in the back-to-back limit. Usinḡ Q ∼ P ⊥ and r k bb /r 2 bb ∼ ∇ k q ⊥ /∇ 2 q ⊥ in Fourier space, we find that the leading contribution of diagrams R 1 × R 2 plus R 1 × R 2 behave parametrically like dσ λ=L This concludes our proof of the power suppression of this diagram relative to the leading power terms. As highlighted in the introduction of this appendix, the odd tensor structure of the hard factor H λ=L,ijk R 1 ×R 2 (P ⊥ ) implies that the contribution computed in this section would vanish at cross-section level, after including the quark-antiquark exchanged diagrams. To get a non-vanishing (but power-suppressed like max(q 2 ⊥ , Q 2 s )/P 2 ⊥ ) contribution at crosssection level, one would need the next term in the small r ⊥ , u ⊥ or u ⊥ expansion. There are several ways to get an even contribution under q ↔q exchange: one of them is to include O(u 2 ⊥ ) corrections from the LO diagram in the complex conjugate. The physical and mathematical arguments presented in this section to demonstrate that diagrams R 1 × R 2 and R 1 × R 2 contribute beyond leading power in q ⊥ /P ⊥ extend to all real diagrams where the gluon interacts at least once with the shockwave -in particular, the diagram R 1 × R 1 . Thus in the back-to-back limit, such diagrams do not contribute to the hard factor but only to the rapidity evolution of the WW gluon TMD.

C Sudakov form factor with two loop running coupling
The Sudakov form factor computed in this paper reads (cf. Eq. (5.36)) with the hard scale µ 2 h = P 2 ⊥ +Q 2 , (C.2) chosen to avoid large negative double logarithms of Q 2 /M 2 qq in the impact factor. The single log coefficients s 0 ands f read s 0 = ln 2(1 + cosh(∆η 12 ))

C.1 Fixed coupling
For fixed coupling α s (µ 2 ) α s (µ 2 h ), we find with the large logarithm in the resummation L ≡ ln and β 1 = . Performing the µ 2 integral, the Sudakov factor simplifies intõ where

D Useful integrals
The following two integrals are useful in computing the LO hard factor for longitudinally polarized photons. They also appear in the LO complex conjugate amplitude in the calculation of the NLO hard factors associated with virtual corrections: (D.1) Differentiating w.r.t. P ⊥ , one gets Additional useful integrals in appendix A.1 are as follows: Differentiating w.r.t. P ⊥ , we find This integral is useful in the calculation of the hard factor for the dressed self-energy: