Dijet photoproduction at low $x$ at next-to-leading order and its back-to-back limit

We compute the cross section for the inclusive photoproduction of a pair of jets at next-to-leading order accuracy in the Color Glass Condensate (CGC) effective theory. The aim is to study the back-to-back limit, to investigate whether transverse momentum dependent (TMD) factorization can be recovered at this order. In particular, we focus on the large Sudakov double logarithms, which are major ingredients of the TMD evolution. Interestingly, the kinematical improvement of the low-$x$ resummation scheme turns out to be a key ingredient in the analysis.


Introduction
A key factor in the success of perturbative Quantum Chromodynamics (pQCD) is the resummation of large logarithms that would otherwise spoil the perturbative expansion. Generally speaking, such logarithms are sensitive to the available phase space for gluon radiation. In one of the most common approaches, known as collinear factorization, the cross section is, up to power corrections, written as a convolution of perturbative hard parts and nonperturbative parton distribution functions (PDFs), and large logarithms in ln(µ 2 /Λ 2 QCD ) are absorbed into the latter with the help of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [1][2][3] evolution equations. A hard scale µ 2 is required to be present in a particle collision for pQCD to be applicable, and is typically provided by the photon virtuality in deep-inelastic scattering (DIS) or by the mass of the produced boson in proton-proton collisions.
Another situation where collinear factorization breaks down is when the process is sensitive to a second scale µ 2 b that is much smaller than the hard scale: µ 2 µ 2 b Λ 2 QCD . Large 'Sudakov' logarithms ln(µ 2 /µ 2 b ) [26] need to be resummed in addition to the collinear ones through the Collins-Soper-Sterman (CSS) evolution equations [27,28] and can be absorbed [29,30] into the parton distributions, which need to be extended to transverse momentum dependent parton distributions functions (TMD PDFs [31]). Classic examples are Z-boson hadroproduction (Drell-Yan) or semi-inclusive deep-inelastic scattering (SIDIS [32]) at low transverse momenta, where the large scale is provided by the boson mass or the photon virtuality, and the small scale by the transverse momentum of the produced boson or hadron.
We should stress, however, that TMD factorization goes beyond Sudakov resummation. Indeed, the quark-and gluon TMD PDFs parameterize the transverse-momentum and spin dependence of the partons inside the proton or nucleus, and moreover come in different types according to the underlying hard process [56]. By contrast, the HEF and BFKL frameworks depend on a single unintegrated gluon distribution and, because of the particular 'nonsense' polarization tensor used, seem to be incompatible with the full structure of the TMD PDFs except at large transverse momenta [57,58].
It turns out that the CGC, which generalizes BFKL, at leading order (LO) in perturbation theory also encompasses TMD factorization. Therefore, one might hope that the CGC provides a unified framework for both TMD factorization and (non)linear low-x evolution. Even at leading order this is not trivial, because a generic CGC cross section involves a complicated intertwining of perturbative coefficients with nonperturbative correlators of semiclassical fields. In the seminal papers [59,60] it was demonstrated that, for a large class of 2 → 2 processes, the CGC and TMD LO calculations do result in the same cross sections, given a proper identification of the correlators of semiclassical gluon fields and gluon TMD PDFs [56,61,62]. This triggered a series of studies, demonstrating the sensitivity to the linearly polarized gluon TMD PDF when masses are included [57,63,64], applying JIMWLK evolution to gluon TMDs [65][66][67], and extending the CGC-TMD correspondence to 2 → 3 processes [68][69][70]. In parallel, the so-called small-x improved transverse momentum dependent (ITMD) factorization framework [58,[71][72][73][74][75][76][77] was developed as a way to combine the applicability of TMD factorization with the resummation of powers (Q s /µ b ) n and (µ b /µ) n , where Q s is the saturation scale, of the CGC.
The aim of this paper is twofold. First, we contribute to the effort to bring CGC calculations to higher perturbative accuracy by calculating the full NLO cross section of inclusive dijet photoproduction, i.e. the process γ +A → dijet +X, using light-cone perturbation theory (LCPT) [78][79][80]. This process could be measured at low photon virtualities in the future Electron-Ion Collider [81], the proposed LHeC [82], or in ultraperipheral leadproton collisions at the LHC. Moreover, our calculation provides an important cross-check of the γ * T + A → dijet + X impact factor recently obtained in [83] using the covariant formulation of the CGC. Second, we want to address the important open question whether the compatibility of the CGC with TMD factorization is preserved beyond leading order.
To do so, we study the limit in which the dijets are back-to-back in the transverse plane, thus creating a scale hierarchy s P 2 ⊥ k 2 ⊥ , where P ⊥ is the typical large transverse momentum of each jet and k ⊥ is their small momentum imbalance. We can reproduce the large Sudakov double logarithms that are essential ingredients in the CSS evolution of the gluon TMD, obtaining the same result as what was predicted in ref. [48,49,84,85]. However, we show that the usual subtraction of low-x logarithms and their absorption into JIMWLK leads to an oversubtraction incompatible with the extraction of the Sudakov logarithms performed in [48,49], and demonstrate that one must rather employ the kinematically-improved JIMWLK equation [86][87][88][89][90][91][92][93][94]. Finally, we observe that, at least at first sight, a class of virtual diagrams which contribute to the finite NLO corrections seem to break factorization. The analysis of these contributions and thus the answer to whether the CGC-TMD correspondence holds at full NLO accuracy is left for future work.
The paper is organized as follows. In section 3, we use light-cone perturbation theory to calculate the loop corrections to γ + A → q +q + X. One set of diagrams: the initialstate loop corrections, were already calculated in ref. [100] and are not recomputed here. In section 3.2, we revisit the calculation of the real NLO corrections, i.e. the process γ + A → q +q + g + X, which was calculated earlier in [69].
As in any NLO calculation, individual diagrams can be plagued by ultraviolet (UV) divergences, while squared diagrams or interferences might exhibit soft and collinear divergences. We demonstrate their cancellation in sections 4, 5, and 8, respectively. Large rapidity logarithms are absorbed in the JIMWLK equation for the LO cross section, as is shown in section 6. The full NLO cross section is then presented in section 9, after which we investigate the back-to-back or 'correlation' limit in section 10.

Leading-order calculation
Throughout this work, we will work in LCPT in the conventions of Bjorken-Kogut-Soper [78,79]. In this picture, the dynamics of the 'projectile', i.e. the photon splitting into the quarkantiquark pair, take place on a much longer timescale than the partonic dynamics of the 'target' proton or nucleus [115,116]. The Color Glass Condensate effective theory then asserts that the target effectively behaves as a localized 'shockwave' of semiclassical gluon fields, which may be described by an external potential built from Wilson lines.
Taking the incoming photon to travel along the + light-cone (LC) direction, colliding head on with the hadronic target which travels along the − LC direction, the differential cross section for the process γA → qqX (see fig. 3.1) is given by: (2.1) The vectors p 1 ≡ (p + 1 , p 1 ) and p 2 ≡ (p + 2 , p 2 ) describe the + and transverse components of the quark and antiquark, respectively, and q + is the photon + momentum. The total + momentum is conserved, as encoded in the delta function. Note that the factor 1/(D − 2) accounts for the averaging over the photon polarization, and that we work for the moment in D dimensions.
In the above formula, the amplitude M is defined as: In the l.h.s. the operatorF , which describes the interaction with the shockwave, acts on the Fock states of the incoming (i) photon and outgoing (f ) quark-antiquark pair. Since these Fock states are asymptotic, their evolution to and from x + = 0: the light-cone time when the scattering takes place, should be calculated up to a given perturbative order: where U is the LC-time evolution operator. At lowest non-trivial order, the photon splits into a quark-antiquark pair before the scattering off the shockwave, while the outgoing quark pair evolves to asymptotic states without perturbative modifications. We thus have from the LCPT Feynman rules [69,79,100]: (2.4) In the above, we introduced the notation PS for the measure of the phase space integrations: where q ≡ (q + , q) and where θ is the Heaviside step function defined as θ(x ≥ 0) = 1 and θ(x < 0) = 0. Combining eqs. (2.3) and (2.4), we can reshuffle the terms in the following form: (2.6) The last term in the above expression encodes the scattering of the 'bare' qq state off the target. In the eikonal approximation, it can be written as: where we introduced the short-hand notation x = µ D−4 d D−2 x, and where denotes a Wilson line in the fundamental representation going from x + = −∞ to x + = +∞ with x − = 0 and transverse position x. We will use the notation W x for Wilson lines in the adjoint representation, and suppress the fundamental color indices. Note that in the eikonal approximation, there is no exchange of spin nor + momentum between the projectile and the target.
Collecting the delta functions in eqs. (2.6) and (2.7), the phase space integration can be written as follows: with the convention q = µ 4−D d D−2 q/(2π) D−2 (note the factor (2π) D−2 which is not present in the integrations over transverse coordinate space). Suppressing the spinor indices, the Dirac structure in eq. (2.6) can be rewritten in function of good spinors [100] with the help of the following intermediary result: which holds irregardless of whether we work with quark or antiquark spinors (since we work with massless quarks). A very useful feature of the above formula is that the good spinors only depend on the + component of the momenta, which means that they are not affected by the shockwave. Note as well that, to arrive at eq. (2.10), we choose the transverse polarization vectors to be i λ = δ iλ . We thus get: where we defined z ≡ p + 1 /q + and σ ij ≡ i 2 [γ i , γ j ]. Note that, in our frame, the photon does not have any transverse momentum.
On the other hand, the energy denominator in (2.6) gives: Putting everything together, we obtain the intermediary result (using the short-hand with: (2.14) Finally, we can perform the integration over p 1 , which gives: Putting everything together, the LO scattering amplitude finally reads: Hence, the amplitude squared becomes: where the brackets denote the average of the semiclassical gluon fields in the target (see also sec. 6.1). The Dirac trace is performed as follows: where we used the completeness relations for good spinors [100]: with the same relation holding for antiquark spinors v G , and with P G = γ − γ + /2 the projector on the good components of the spinor field.
Since P G commutes with all transverse gamma matrices and since P G P G = P G , we get: (2.21) We finally obtain, using the identities (A.5) and (A.7): Our final result for the LO amplitude squared is then: which leads to the following cross section 1 in D = 4 dimensions: (2.24) In the above two expressions, we introduced the following compact notations for the dipole and quadrupole: (2.25) 3 Next-to-leading order diagrams The calculation of the next-to-leading order amplitudes follows largely the same method as the leading-order case. We will therefore not reproduce all the intermediate steps but Real NLO Table 1. Overview of the divergencies and large logarithms encountered in our calculation, for the different real NLO contributions to the cross section. In the column 'cross', all cross terms are meant between the gluon emissions before and after the shockwave, from the quark or from the antiquark, i.e. the interference between M QFS , M QFS , M QSW , M QSW . Terms involving a real instantaneous gluon emission are strictly finite and do not contribute to the Sudakov logarithms at our accuracy. In this work, we do not attempt to analyze the Sudakov double logarithms beyond the large-N c limit. Moreover, in our regularization scheme, it is not always possible to unambiguously distinguish soft from rapidity divergences. The question mark indicates this is the case for |M QSF | 2 and |M QSF | 2 . The only certainty is that these two diagrams combine with FSIR (table 2) into a contribution to the cross section that has rapidity divergences only, see also section 8.
rather give the result for each Feynman diagram, depicted in figures 3.1 and 3.3. Of course, all diagrams have a counterpart in which the quark and antiquark are reversed. We will denote these contributions with an overline: for example, GESW is the graph in which the gluon is radiated by the quark and, after scattering off the shockwave, absorbed by the antiquark. These 'q ↔q conjugate' amplitudes can be obtained in a straightforward way from their counterparts applying the following steps: 1. Interchange all the indices 1 and 2 except in the Wilson lines and in the spinors u G (p + 1 ) and v G (p + 2 ).
2. Take the complex conjugate of the part of the Dirac structures sandwiched between the spinorsū G (p + 1 ) and v G (p + 2 ).
3. In LCPT, the vertex for the emission or absorption of a gluon from the antiquark has an overall minus sign. Add it to the diagrams ISW, QSW and QSF.
4. Calculate the Wilson line structure separately, there is no simple rule here.
In tables 1 and 2, we list all NLO contributions to the cross section and their possible contribution to ultraviolet or infrared divergencies. Their regularization and either cancellation or renormalization will be the subject of sections 4 to 8. The amplitude for the instantaneous production of a quark, antiquark, and gluon from the photon, where the gluon crosses the  Virtual NLO   SESW, sub GESW GEFS + IFS SESW, UV IS FSIR FSUV  ultraviolet     JIMWLK     ?  ?   ?  collinear   soft   ?  ? ? ?  shockwave and is absorbed by the outgoing quark ( fig. 3.1), is found to be:  with

Virtual corrections ISW: Instantaneous gluon traversing SW
and where we introduced the generalization of the Coulomb field (see [69]) to D − 4 dimensions, defined as: ig e e f g 2 To arrive at the above expression, we made use of the intermediary result: where A is the modified Weizsäcker-Williams field 2 defined as: (3.6) Luckily, it turns out that we will only need to evaluate A in finite contributions to the cross section, where we can work in D = 4 dimensions for which its expression reduces to: (3.7) The spinor structure in (3.4) has the following form: where we defined: which is the most singular part of the Dirac structure, scaling like 1/(k + 3 ) 2 .
SESW: Self-energy correction traversing the SW For the self-energy corrections crossing the shockwave (SESW, fig. 3.1), we obtain: (3.10) The above amplitude contains an ultraviolet divergence, which comes from the limit where x 3 → x 1 . Following ref. [100], we construct a counterterm M SESW,UV by taking the x 3 → x 1 limit in M SESW except in the singular part, and by subtracting an infrared (IR) divergent contribution ∝ Aη(x 13 )Aη(x 23 ) as follows: where we used that, even in D dimensions: By construction, the counterterm M SESW,UV has the same UV pole as M SESW while not possessing any other divergences. To show this is the case, let us explicitly evaluate the x 3 integration in eq. (3.11). We have that: since scaleless integrals disappear in dimensional regularization. This is equivalent to the statement that the above integral contains two divergences, an UV (x 3 → x 1 ) and an IR (x 3 → ∞) one, which cancel each other. The second integration in eq. (3.11) reads: 14) which is divergent in the infrared. Combining eqs. (3.13) and (3.14), the IR pole of the latter cancels the one hidden in the former, and the overall divergence of M SESW,UV can be interpreted as an UV one. Since M SESW,sub ≡ M SESW − M SESW,UV is now free from UV divergences, we can evaluate it in D = 4 dimensions: with: In the above expression, P ⊥ is a transverse vector defined as: while the loop momentum K is related to the virtual gluon transverse momentum through IFS: Instantaneous gluon exchange in the final state In the case of an instantaneous qq → qq final-state interaction mediated by an instantaneous gluon ( fig. 3.1) we obtain: where the loop momentum is again defined as p Combining diagrams GEFS and IFS The amplitudes (3.16) and (3.19) both exhibit an unphysical (1/k + 3 ) 2 power divergence, which cancels when summing them. In diagram GEFS, this divergence stems from the Dirac structure Diracλη η q→q,(ii) . The subamplitude corresponding to this part of the Dirac structure, which we denote by M GEFS,(ii) , reads: 20) and is clearly divergent ∝ 1/(k + 3 ) 2 in the k + 3 → 0 limit. The situation changes when the above (sub)amplitude is summed with the IFS amplitude (3.19): and we are left with a logarithmic k + → 0 divergence, which will contribute to JIMWLK.
Final state with gluon in s-channel Finally, in fig. 3.2, two virtual diagrams are depicted in which, after the scattering off the shockwave, the quark-antiquark pair annihilates and a new pair is created. This process either takes place through an s-channel gluon or as an instantaneous qq → qq vertex mediated by a fictitious instantaneous gluon. When multiplied with the conjugate leading-order amplitude, these two diagrams contribute to the cross section with a coupling: In the above, the summation over quark flavors takes place for the two quark lines separately, since at the level of the cross section the two fermion loops are distinct and only connected by a gluon. This is a unique feature of the two diagrams under consideration; all other contributions to the NLO dijet photoproduction cross section have a coupling f e 2 f g 2 e g 2 s . In particular, since in this calculation we only consider the three lightest quarks: hence these diagrams can be disregarded.
with the Dirac structure: which can be written as: (3.27) RI: Real gluon created instantaneously before the shockwave Finally, the real gluon can be radiated instantaneously from the photon in addition to the quark-antiquark pair ( fig. 3.3), yielding the amplitude: where (3.29)

UV safety
The diagrams SESW and SESW, i.e. the quark-and antiquark self-energy corrections where the gluon crosses the shockwave, are UV divergent and need to be regularized by adding a counterterm, as we demonstrated in eqs. (3.10)-(3.15). These are the only diagrams which we calculated that exhibit a UV singularity. There are, however, virtual contributions that we did not explicitly compute but whose UV poles will cancel with the one in our counterterm, resulting in a UV-finite total NLO amplitude. This is what we will now demonstrate. First, comparing expression (3.11) for the UV counterterm with the LO amplitude (2.17), we can write: where, in the last line, for future convenience we introduced a slight abuse of notation since the factorization of the term V SESW,UV is really on the integrand level. This term reads: To arrive at the above expression, we made use of eqs. (3.13) and (3.14) and introduced the following phase space integral over the gluon + momentum: where k + min is a regulator which will be further specified in subsection (6.1). Writing D = 4 − 2 UV , one can calculate further to obtain: and: such that: Likewise, the contribution from the antiquark self-energy loop scattering the shockwave gives: hence in total: In ref. [100] (see eq. 144), it was demonstrated that the loop corrections to the initial state, which we did not explicitly calculate in this work, can be cast into a similar form factor V IS as the one that appears in (4.1): Adding eqs. (4.8) and (4.9) gives: (4.10) Clearly, adding the initial-state loop corrections is not yet sufficient to cancel the UV poles from the SESW counterterms, and one needs another UV-divergent contribution. This final ingredient is provided by the self-energy corrections to the quark and antiquark in the final state. We did not compute these diagrams explicitly, since they are simply zero. This is because in dimensional regularization the UV and IR contributions to a scaleless integral -such as self-energy loops on an asymptotic massless state-exactly cancel (cfr. eq. (3.13)). This property can, however, be exploited to construct the missing piece for the UV cancellation in (4.10), writing the total contribution of the asymptotic (anti)quark leg corrections as V FS = V FSUV + V FSIR = 0 with: where we added subscripts to distinguish between positive infinitesimal numbers parameterizing the UV resp. IR pole. Therefore, the sum of the leading-order diagram (2.6), the initial state loop corrections (4.9), the UV counterterms from the SESW and SESW diagrams (4.8), and the UV-divergent part of the final-state corrections (4.11), is UV-finite: (4.12) After this procedure, what is left from the cancellation of UV divergencies is the finite term above, as well as a new contribution to the virtual diagrams which stems from the IR-divergent part of the asymptotic leg corrections (4.11): which ultimately will cancel with the IR poles from the real NLO corrections (see section 8).

Soft safety in gluon exchange and interferences
A second type of infinities that appear in our calculation are soft divergences, which show up when the gluon momentum k 3 = (k + 3 , k 3 ) → 0. They are distinct from the so-called rapidity divergences associated with k + 3 → 0 but where k 3 stays finite. The latter will be regulated with a cutoff k + min and the remaining logarithms resummed by the JIMWLK equation (see section 6).

Virtual contributions
When computing the virtual diagrams in section 3.1, we have performed the integration over the transverse loop momentum k 3 . Therefore, the information on possible soft singularities has been lost. To investigate whether a diagram contains a soft divergence, one needs to rescale the gluon transverse momentum with its + momentum: k 3 →k 3 = (k + 3 /p + 1 )k 3 and then take the k + 3 → 0 limit. The only two virtual diagrams with a soft singularity are the ones with a gluon exchange in the final state: GEFS and IFS (and their q ↔q counterparts).
Let us start with the GEFS diagram and perform a rescaling K = χ /z with Performing a first-order Taylor expansion in denominator, which tends to zero when χ → 0, one obtains in the limit: The only other source of powers 1/k + 3 in M GEFS (3.16) is the Dirac structure, which needs to scale as 1/(k + 3 ) 2 in order to combine with the above integral to give a singularity. We can therefore pick up only the simple contribution Diracλη η q→q,(ii) to the full Dirac structure (3.9), which gives: In the above expression, Diracλ LO is the q ↔q conjugate of the leading-order Dirac structure eq. (2.14)), obtained by interchanging 1 ↔ 2 and taking the complex conjugate of the structure between the spinors (cfr. the algorithm in the beginning of sec. 3): Combined with the delta function from the Dirac structure, the soft limit of the integral J in (5.2) is thus reduced to two integrals which both can be analytically solved in dimensional regularization: Let us start with the integral I 2 . Applying the 'real' resp. 'complex' Schwinger trick to the first and the second denominator: we obtain Evaluating the integral over α, one finds: Using the integral representation of the gamma function: we end up with: Writing −i = e i(− π 2 +2nπ) and D = 4 − 2 IR : Likewise, applying (5.6) to the (finite) integral I 1 yields: The next step is to rewrite the D − 2 dimensional -integration as follows: (5.13) The I 1 integral can then be evaluated starting with the integration over cos θ, followed by the β integral: (5.14) -20 - We finally obtain: The soft limit of diagram IFS can be extracted in a similar way. Rescaling The IR pole will be canceled with certain real NLO corrections, as will be shown in the next subsection. This cancellation takes place on the level of the amplitude squared, obtained by multiplying with M † LO . Adding as well the q ↔q diagrams and the complex conjugate (c.c.), we obtain: where we made use of the spinor trace (2.22) with D = 4, as well as:

Real contributions
The IR pole found in the previous subsection stems from a soft virtual gluon exchange in the final state. It will cancel with real contributions that have the same topology, i.e. a soft final-state gluon radiated from the quark in the amplitude and from the antiquark in the complex conjugate amplitude, or vice versa. The corresponding contribution to the cross section is: and taking the p + 3 → 0 limit, the above equation becomes: where we wrote dp + 3 /p + 3 → dξ/ξ and extracted the leading-power of the Dirac structure: which leads to: The integral over transverse momentum can be cast into the following form (where P ⊥ is again the momentum vector defined in (3.18)): (5.24) The first integral disappears in dimensional regularization because it is scaleless. The second one can be computed by applying the real Schwinger trick (5.6) twice: where we wrote −1 = e i(π+2nπ) . Combining the above result with eq. (5.6)) and adding the complex conjugate, we finally obtain: The above result is exactly the opposite of the soft limit of the virtual contributions, eq. (9.3). Therefore, the total cross section is free from soft divergences from contributions with final state gluon exchange topology (including interferences of real final state emission from the quark or the antiquark). Additional soft divergences, appearing together with collinear divergences, will be discussed in section 8.

Kinematics
So far, among the virtual NLO corrections to the dijet cross section that we have calculated, many have a logarithmically divergent integral over the + momentum k + 3 of the virtual gluon stemming from the k + 3 → 0 regime. In some cases, like the SESW, UV contribution (3.11) or the IS contribution (4.9), in which the transverse loop integration can be performed explicitly, one could have used dimensional regularization to deal with these divergences. But in other cases, like the GESW contribution (3.4) or the SESW, sub contribution (3.15), the integration over the transverse position of the gluon x 3 cannot be performed explicitly due to the presence of a Wilson line at x 3 . In such cases, dimensional regularization cannot handle the k + 3 → 0. For this reason, we introduce the lower cutoff k + min to regularize all k + 3 loop integrals. Similarly, one encounters divergences in the real NLO corrections at the cross section level from the p + 3 → 0 regime in the integration over the real gluon + momentum p + 3 (either inside or outside the measured jets, as we will see in sections 7 and 8), which we regularize with the same cutoff k + min . Part of these k + 3 → 0 or p + 3 → 0 divergences are genuine soft divergences which cancel at the dijet cross section level, as we have seen in the previous section 5. However, others are rapidity divergences which survive in the form of single logs of k + min after regularization. These large logs of k + min are high-energy leading logs, which can be extracted from the NLO correction to the cross section and resummed into the LO term, via the JIMWLK evolution of the target-averaged color-or Wilson-line operator, as we will now explain.
In the LO cross section (2.24), the target-averaged color operator is unevolved and should not yet include high-energy logarithms. We can, therefore, use the notation: In the simplest scheme for the JIMWLK evolution, which we will use in most of the present study, JIMWLK is viewed as an evolution equation along the k + axis in logarithmic scale.
In this scheme, one defines as the same target-averaged operator but now including the resummation of high-energy leading logs associated with gluons with light-cone momentum k + between the cutoff k + min and the factorization scale k + f , in the notations of ref. [91]. The evolution with the factorization scale k + f , or equivalently with Y + f , is given by the JIMWLK equation for the LO operator A more explicit version of this equation, with the action of the JIMWLK Hamiltonian H JIMWLK fully worked out, can be found in ref. [117]. Integrating eq. (6.3), one finds The JIMWLK Hamiltonian is of order α s . Hence, the dependence of a target-averaged operator on Y + is an effect suppressed by one extra power of α s in fixed-order perturbation theory. Therefore, expanding eq. (6.4) in powers of α s we can write with the scale unspecified for the operator in the second line, since it is not under control at this perturbative order. Hence, inserting eq. (6.5) into the LO cross section (2.24), one substitutes the unevolved target-averaged operator with its evolved (up to the factorization scale k + f ) version, generating an extra NLO term which involves the JIMWLK Hamiltonian. This new NLO term will subtract the logarithmic dependence on the cutoff k + min found in the NLO cross section due to rapidity divergences at k + 3 → 0 or p + 3 → 0. Writing formally the NLO correction to the dijet cross section found from the fixedorder calculation as to separate the gluon + momentum integral from the rest of the cross section (with other Heaviside or Dirac delta functions constraining p + 3 included inσ NLO ), one has: By construction, the first term in eq. (6.7), extracted from the total NLO correction, identically cancels the second term from eq. (6.5) after substituting the left hand side of eq. (6.5) into the LO cross section (2.24). Then, the statement that rapidity divergences are subtracted and resummed thanks to the JIMWLK evolution is equivalent to saying that in the second term of eq. (6.7), the cutoff k + min can be dropped, thanks to cancelations happening at low p + 3 between the terms in the square bracket. In the rest of this section, we will check this statement, by studying the p + 3 → 0 (or k + 3 → 0) limit for each of the NLO contributions to the cross section.
Finally, we should discuss appropriate values for k + f and k + min (and thus for Y + f ) in order to resum high-energy logarithms via JIMWLK evolution. The factorization scale k + f should be of the order of the + momenta of the measured jets (or at most q + ). The cutoff k + min represents the typical + momentum scale set by the valence (and other large-x) partons inside the target, before evolution. Modeling the target before low-x evolution as a collection of partons carrying a fraction of at least x 0 of the target momentum p − A and with a typical transverse mass Q 0 , one has: Moreover, we have chosen a frame in which the photon momentum q µ = (q + , 0, 0) lies entirely in the light-cone + direction, and the target nucleus where M A is the target mass. Then, the total energy squared of the collision is at high energy. Hence, one can write the cutoff as x 0 s , (6.10) and the range for JIMWLK evolution as For this reason, Y + f is considered to be a high-energy logarithm. x 0 can be taken to be 0.01, or at most 0.1. Q 0 should be a scale around the transition between perturbative and non-perturbative QCD, or should be related to the initial saturation scale Q s,0 in the case of a large enough nucleus. However, in practice, Q 2 0 /x 0 can be treated as a parameter in a BK/JIMWLK global fit, together with the shape of the initial condition for the evolution.
The scheme chosen for the JIMWLK resummation, based on an evolution strictly along the p + axis, is particularly simple to handle. However, this scheme is not unique and, in fact, neither is it optimal as we will discuss in section 10.3, in particular for the study of Sudakov logarithms.
In the rest of this section, we start by considering the virtual NLO amplitudes and studying their k + 3 → 0 limit. After that, we bring them to the level of the cross section by multiplying them with the complex conjugate of the LO amplitude M † LO , constructing the total virtual contribution to JIMWLK. For the real NLO amplitudes the procedure is similar, as it turns out to be easiest to take the p + 3 → 0 limit at the amplitude level. Interestingly, we find that the thus obtained 'virtual' and 'real' contributions to JIMWLK are separately free of subleading-N c terms. In the end, we demonstrate how in the k + 3 , p + 3 → 0 limit, the cross section corresponds to the JIMWLK evolution equations applied to the Wilson-line structure Q 122 1 − s 12 − s 2 1 + 1 (6.12) of the leading-order result, which confirms the resummation of high-energy logs by JIMWLK into the LO term, as presented in this section.

Virtual diagrams
GEFS+IFS In the k + 3 → 0 limit, the subamplitude M GEFS,(ii)+IFS becomes: There is another contribution to JIMWLK due to the Diracλη η q→q,(i) spinor structure in the amplitude M GEFS , which yields: (6.14) Subamplitudes eq. (6.13)) and (6.14) nicely combine into: Finally, with the help of definition (2.16) of the Weizsäcker-Williams fields, it is easy to show that: Multiplying with the complex conjugate of the leading-order amplitude, we finally obtain: GESW Since the modified Weizsäcker-Williams structure is finite in the limit k + 3 → 0: the only contribution to JIMWLK from this diagram comes from the Dirac q→q(ii) term: After multiplying with M † LO , making use of eq. (5.19)), one obtains: It is easy to see that the q ↔q counterpart of this diagram will give the contribution: SESW Taking the k + 3 → 0 limit of (3.15) is trivial and yields, after multiplying with M † LO : Likewise, we get for the diagram with a gluon loop on the antiquark: FSIR The last set of virtual diagrams that exhibit a rapidity divergency and hence contribute to JIMWLK are the IR parts of the self-energy corrections to the asymptotic (anti)quark, eq. (4.13)): where, from eq. (4.11)): Total virtual contribution to JIMWLK Collecting all the above virtual contributions to JIMWLK, we finally obtain: An interesting feature of the above formula is that all subleading-N c contributions have cancelled.

Real diagrams
Taking the p + 3 → 0 limit of the real gluon emission amplitudes M QFS and M QSW , eq. (3.27) resp. (3.24), is very straightforward due to the simple leading-power behavior (5.22) of the Dirac structure: Care should be taken with the amplitudes M QFS and M QSW , which as we remarked in sec. 3, receive an additional minus sign from the LCPT Feynman rules due to theq →qg vertex: lim (6.31) With the above expressions at hand, constructing the real part of the JIMWLK equation is a trivial task, which yields:

(6.33)
The above expression is the JIMWLK equation for M LO 2 , eq. (2.23)), consisting in the evolution of the quadrupole Q 122 1 (eq. 4 in Ref. [117]) and, in the last line, of the dipole s 12 and its complex conjugate. We have therefore proven what we asserted in eq. (6.7)), namely that the part of the cross section dσ enhanced by large logarithms ln k + f /k + min takes the formĤ JIMWLK dσ LO .
The structures A , B , K 1 and K 2 in eq. (6.33)) are defined according to the notation in [117], and read: , (6.35) , (6.36) In the previous sections, we were always concerned with partonic scattering amplitudes, hiding any hadronic physics inside the nonperturbative CGC averages over the Wilson lines. The aim of this work, however, is to compute the NLO dijet photoproduction cross section, not the diquark one, which implies that at a certain point we need to quantify what we mean by a jet. At leading order, transforming the γA → qqX cross section to the γA → dijetX one is trivial and done by simply identifying the quark and antiquark partonic momenta p 1 and p 2 with the jet momenta p j1 and p j2 . This same trivial identification can be done for the virtual NLO contributions to the cross section. It is implicitly assumed that the outgoing momenta are sufficiently separated to prevent the quark and antiquark being grouped within the same jet.
For the real NLO corrections, the situation is different. In this case, two steps are needed to go from the tree-level γA → qqgX partonic cross section to the real NLO correction to the dijet cross section. First, a jet algorithm needs to be applied in order to determine in what part of the phase space the three partons are considered to form three separate jets, and in what part of the phase space two of the partons are clustered into a single jet. Second, in the case of the three jet configuration, in order to obtain the corresponding contribution to the inclusive dijet cross section at NLO, two jets should be identified with the measured jets with momenta p j1 and p j2 while the third jet should be integrated over.
The jet definition/algorithm that we use in this NLO calculation is as follows: if two partons i and j in the final state are such that then they are considered to be part of the same jet with momentum (p + i + p + j , p i + p j ). Otherwise, the partons i and j are considered to form separate jets. This definition depends on the jet radius parameter R, satisfying 0 < R < 1.
In the present study, for simplicity and in order to be able to perform analytic calculations as far as possible, we are using our jet definition in the narrow-jet limit, meaning that formally R → 0. In practice, this means that terms in log(R) or terms independent of R are kept in the cross section, whereas terms suppressed as positive powers of R are neglected. In diagrams without a collinear divergence, one can identify each of the three partons with a separate jet, because the contribution from the merging of two partons into one jet takes place only in a parametrically small part of the phase space, suppressed by powers of R. Among the real NLO corrections, only |M QFS | 2 has a quark-gluon collinear divergence and only |M QFS | 2 has an antiquark-gluon collinear divergence. In the narrow-jet limit, it is then sufficient to distinguish when the quark and the gluon (resp. the antiquark and gluon) are clustered into a single jet or not by the jet algorithm in the |M QFS | 2 contribution (resp. in the |M QFS | 2 contribution).
As a remark, when the left hand side of eq. (7.1) is much smaller than 1, one has the equivalence where ∆y ij and ∆φ ij are the difference in rapidity and in azimuthal angle between the two partons. For this reason, our jet algorithm and the Cambridge/Aachen algorithm (see refs. [118,119]) become equivalent in the narrow-jet limit R → 0. Moreover, in the case of three-jet configurations, we always consider the quark jet to be the first measured jet, the antiquark jet to be the second measured jet, and the gluon jet to be the jet which is unmeasured and integrated over. Hence, the results that we present correspond specifically to the flavor-tagged inclusive quark-antiquark dijet cross section at NLO. Other contributions to the full inclusive dijet cross section at NLO, for example with a gluon jet and a quark jet, do not involve any conceptual difficulty (in particular, they do not contain any divergence whatsoever), and could be obtained in a straightforward way from our intermediate results on the quark-antiquark-gluon partonic cross section.
From the |M QFS | 2 contribution to the qqg partonic cross section applying our jet definition in the narrow-jet limit, we thus obtain two contributions to the dijet cross section. The first one, in which each of the three partons is forming a jet and in which the gluon jet is integrated over, reads dσ |QFS| 2 ; out dijet dp + j1 d D−2 p j1 dp + j2 d D−2 p j2 We have excluded the region in which the quark and the gluon are merged into one jet thanks to θ in ( p i , p j ), which is defined as and enforces the condition (7.1). The second contribution, in which the quark and the gluon are combined into one jet of momentum p j1 = p 1 + p 3 whereas the antiquark forms a jet of momentum p j2 = p 2 , reads dσ |QFS| 2 ; in dijet dp + j1 d D−2 p j1 dp + j2 d D−2 p j2 For the |M QFS | 2 contribution, which is the square of the gluon emission by the antiquark in the final state, one obtains two contributions to the NLO dijet cross section in a similar way. Their expressions are the same as eqs. (7.4) and (7.6), up to the exchange of the role of the quark and the antiquark, meaning in particular p 1 ↔ p 2 and p j1 ↔ p j2 .
Finally, for all the contributions to the qqg partonic cross section other than |M QFS | 2 and |M QFS | 2 , which are collinear safe, we go from the partonic to the dijet cross section by integrating over the gluon momentum and identifying the produced quark and antiquark with the measured jets, as dσ NLO real; coll. safe diags. dijet dp + j1 d D−2 p j1 dp + j2 d D−2 p j2 = PS( p 3 ) dσ coll. safe diags.

Collinear and soft safety in final state fragmentation
As already mentioned in the previous section, the real contributions |QFS| 2 and |QFS| 2 to the dijet cross section have collinear divergences. Such divergences come from the transverse integration, and we are handling them with dimensional regularization. In addition, these diagrams also lead to soft divergences at the dijet cross section level. In our regularization scheme, using dimensional regularization for transverse integrals and the cutoff k + min for the p + 3 integral, it is difficult to distinguish unambiguously between genuine soft divergences and rapidity divergences, since both arise from the p + 3 → 0 regime. In this section, calculating the real contributions from |QFS| 2 to the dijet cross section, we will show that most divergences cancel in the sum of the real corrections |QFS| 2 and |QFS| 2 on the one hand, and of the virtual correction FSIR both in the amplitude and the complex conjugate amplitude on the other hand. The only leftover divergence in the resulting sum is the rapidity divergence associated with the JIMWLK evolution, already discussed in section 6. This is an evidence for the cancelation of both collinear and soft divergences between these contributions. In order to calculate the real correction |QFS| 2 to the dijet cross section, following the discussion in section 7, we first split it into the in contribution (7.6), in which the quark and the gluon belong to the same jet, and the out contribution (7.4), in which they form separate jets.
The amplitude for the diagram QFS is given in eq. (3.26). Squaring it and summing over the colors, helicities and the photon polarization, one finds

Contribution |QFS| 2 ; in
When the gluon and the quark are merged into a single jet by the jet algorithm, we have by definition p 1 + p 3 = p j1 . Introducing the transverse momentum of the gluon relative to the one of jet j1 one can rewrite the denominator of eq. (8.1) as Then, comparing with the LO squared amplitude (2.23), eq. (8.1) becomes: Moreover, with the notation (8.2), Inserting eqs. (8.4) and (8.5) into eq. (7.6), one finds the inside jet radiation contribution from |QFS| 2 to factorize as In eq. (8.7), the transverse integral is straightforward to calculate in dimensional regularization, and yields Note that we perform the D → 4 expansion before taking the integration over p + 3 , since in our scheme, only the transverse integral (and thus the UV and collinear divergences) is regulated with dimensional regularization, whereas the p + 3 -integral is regulated with the lower cutoff k + min . Thus, eq. (8.7) becomes: As explained in section 6.1, k + min plays a double role in our calculation. On the one hand it is used as a regulator for the integrals in p + 3 or k + 3 . On the other hand, it is used to specify the physical scale set by the target at which the low-x evolution is starting, or equivalently to encode the dependence on the total energy √ s of the collision. The JIMWLK evolution resums single high-energy logarithms, written as logarithms of k + min , and unlike the jet radius parameter R does not depend on details of the process. Hence, in the result (8.9), the term in ln 2 k + min and the term in ln R ln k + min can definitely not be subtracted and resummed by JIMWLK. Instead, these two terms in eq. (8.9) should be understood as manifestations of the standard soft-collinear double logarithmic divergence for final state gluon radiation in our hybrid regularization scheme, having nothing to do with rapidity divergences and low-x evolution. In the remainder of this section, we will see how these ln 2 k + min and ln R ln k + min terms as well as the collinear 1/ coll pole are canceled by other contributions, leaving only the expected single high-energy logarithm to be subtracted and resummed by JIMWLK, following section 6.
For the diagram |QFS| 2 , the contribution from the regime in which the antiquark and the gluon are merged into the same jet can be calculated in the same way, leading to (8.10)

Contribution |QFS| 2 ; out
Let us now consider the other contribution to the dijet cross section from the |QFS| 2 diagram, in which the quark and the gluon are forming separate jets according to our jet definition. In that case, the quark and antiquark momenta are identified with the momenta of the measured jets, as p 1 = p j1 and p 2 = p j2 . Using again the notation (8.2), the expression (8.3) for the denominator in eq. (8.1) stays valid. By contrast, the x 11 dependent phase factor in eq. (8.1) now becomes: With all this, the squared amplitude |M QFS 2 from eq. (8.1) can be written in the case of three separate jets as with M LO 2 written as in eq. (2.23). Note that eq. (8.12) is, once again, a slight abuse of notation, since the factorization actually happens happens at the integrand level and the integrations over x 1 and x 1 are hidden inside |M LO | 2 . Moreover, from the definition (7.5), one has now In the narrow jet R → 0 limit, this theta function changes values from 0 to 1 at a parametrically small value of |P 3 |, which scales as R. Hence, in that limit, P 3 is negligible compared to p j1 , so that Inserting the expressions (8.12) and (8.14) into eq. (7.4), we can write the contribution to the dijet cross section from |QFS| 2 with three separate jets as .
In eq. (8.15), the integral over P 3 is actually finite, since the phase prevents UV divergences to occur, and the theta function is cutting off the collinear regime, preventing the gluon to belong to the quark jet. Hence, dimensional regularization is actually not necessary here and one can take D = 4 − 2 → 4. The transverse integral in eq. (8.15) can be calculated as (8.16) in the narrow-jet limit R → 0. In eq. (8.16), J 0 is the Bessel function of the first kind, and c 0 ≡ 2e −γ E . At this stage, eq. (8.15) becomes dσ |QFS| 2 ; out dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 In this expression, we could simply use the delta function in order to perform the p + In order to obtain a well behaved result and to extract the sensitivity of the expression (8.17) on the cutoff k + min , let us split it into several pieces. First, let us take the integrand at p + 3 → 0, including in the delta function. The corresponding piece of eq. (8.17) can be written as dσ |QFS| 2 ; out; soft dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 introducing p + j1 as an upper bound in this soft contribution. The first two terms in eq. (8.19): in ln 2 k + min and ln R ln k + min , precisely cancel the ones found in from the in-jet radiation contribution (8.9). As discussed earlier, these terms cannot be associated with high-energy evolution but are, instead, manifestations of soft-collinear divergences in our regularization scheme. Their cancellation is a crucial requirement for the consistency of our hybrid regularization scheme in which we combine transverse dimensional regularization with a lower cutoff for p + 3 . Subtracting the term (8.18) from (8.17) is sufficient to remove any divergence from the p + 3 integral, or more precisely any logarithmic sensitivity on k + min , so that one can remove the cutoff k + min in the leftover piece. Nevertheless, the finite p + 3 integral can still produce a potentially large logarithm. The reason for this is that the p + 3 → 0 approximation of the integrand of eq. (8.17) in order to obtain (8.18) is not always a good approximation for small but finite p + 3 p + j1 , p + j2 . Indeed, in this regime, the smallness of p + 3 compared to p + j1 can be compensated for very large values of |p j1 · x 11 |, making the phase factor in eq. (8.17) non trivial even at small p + 3 . In order to have control on this issue, let us extract as well from eq. (8.17) the contribution dσ |QFS| 2 ; out; phase dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 where V out; phase The above contribution will be further studied in section 10, in the back-to-back jets limit.  (8.22) since in this contribution, the integration over p + 3 cannot produce a divergent result even without cutoff and cannot produce further large logarithms. It can be written explicitly as dσ |QFS| 2 ; out; reg dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 For the diagram |QFS| 2 , the contribution from the regime in which the antiquark and the gluon form separate jet can be split in the same way into three contributions, and one obtains V out; soft and dσ |QFS| 2 ; out; reg dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 = dσ |QFS| 2 ; out; reg dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 (1 ↔ 2) . (8.26)

Cancellation of collinear and soft divergences
As already discussed, soft and collinear divergences manifest themselves in various ways when |QFS| 2 is calculated using transverse dimensional regularization and a cutoff k + min . In addition to the 1/ collinear pole in the inside jet radiation contribution (8.9), both eqs. (8.9) and (8.19) contain terms of soft-collinear origin, namely in ln 2 k + min and in ln R ln k + min . It is then difficult to disentangle single logarithms ln k + min associated with either soft or rapidity divergences.
Adding eqs. (8.9) and (8.19) together, the terms in ln 2 k + min and in ln R ln k + min cancel, and one obtains The other two contributions (8.21) and (9.26) from |QFS| 2 do not contain any divergence. Similarly, for the diagram |QFS| 2 , one finds In the NLO virtual corrections to DIS dijets, the only contribution containing a collinear divergence is FSIR, see eq. (4.11), which factorizes as well from the LO cross section at integrand level. That contribution correspond to the IR part of the self energy diagrams for the quark and for the antiquark in the final state. With the identification of the momenta of partons and jets, it can be written as Eq. (8.29) is corresponds to a virtual correction in the amplitude. We should also consider that virtual correction in the complex conjugate amplitude, which is The only difference between the two is the dependence on the transverse coordinates In this expression, all soft and collinear divergences have canceled [120]. Only single logarithmic dependence on k + min remains, which now correspond to high-energy evolution. This can be shown as follows. Introducing the factorization scale k + f , one can isolate the dependence on the cutoff k + min as (8.32) The first term in eq. (8.32) amounts to a multiplication of the LO cross section by This corresponds indeed to the part of the JIMWLK evolution proportional to α s C F times the LO cross section. Hence, in our resummation scheme for high-energy leading logs, the first term in eq. (8.32) is the part which is subtracted and resummed by the JIMWLK evolution.

Inclusive dijet cross section
In this section, we present the full differential NLO cross section for γA → dijetX in the CGC. It is given by the sum: where it is understood that the rapidity divergences are absorbed into JIMWLK according to eq. (6.7). The cross section is written in terms of different contributions which are each separately soft-and collinear safe. We will shortly discuss each of the terms, and provide their explicit expressions in the next subsections. The first part of the cross section is the leading-order one, given in eq. (2.24)), where the outgoing partons are identified with jets: p 1,2 → p j1,j2 .
Second, following the discussion in sections 7 and 8, the real NLO corrections due to final-state gluon emission contain soft-collinear singularities. These singularities are cured when applying the jet algorithm and summing the contributions due to gluon emission inside the jet, soft gluon emission outside the jet, and the IR part of the virtual self-energy corrections to the final state: The explicit expression for V jet is given in eq (8.32). Third, as explained in section 5, the contribution to the cross section due to diagram GEFS, (ii) + IFS and its q ↔q counterpart contains a soft divergence. This singularity cancels with the one found in the interference between real final-state emission from the quark in the amplitude, and from the antiquark in the complex conjugate amplitude or vice versa. Summing both virtual and real contributions, one ends up with the following finite contribution to the cross section which we might call softGE for soft final-state gluon exchange: Fourth, there are terms due to the final-state gluon radiation outside the jet, that in certain kinematics will become enhanced by a large Sudakov double logarithm (see sections 8 and 10): The explicit formulae for V out; phase The final contributions to eq. (9.1) are due to purely finite virtual and real NLO corrections: and dσ finite real dp + j1 d 2 p j1 dp + j2 d 2 p j2 = dσ |QFS| 2 ;out;reg dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 + dσ |QFS| 2 ;out;reg dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 with the Dirac trace, calculated with the help of identity (A.8): GESW We find for the contributions due to a gluon being exchanged between the antiquark and quark, traveling through the shockwave: where the Dirac trace is easily calculated applying identity (A.8), with the following result: The q ↔q conjugate, where the gluon is emitted from the quark and, after interaction with the shockwave, absorbed by the antiquark, is given by: 16) and the Dirac trace yields: ISW The following contributions are due to the instantaneous splitting of the photon into a gluon, quark and antiquark. After traveling through the shockwave, the gluon is absorbed by the quark: where we used the following result for the Dirac trace: with the Dirac trace: Tr Dirac λ † LO Diracλη η q→q,(i) = 32 The q ↔q conjugate contributions to the cross section can be simply obtained from the above result by exchanging the quark and antiquark indices: IS+UV+FSUV Finally, the contributions due to the initial-state corrections, ultraviolet counterterms, and ultraviolet part of the self-energy corrections to the asymptotic final states read:

Real terms
Final-state gluon radiation outside jet The two terms in the first line of eq. (9.6) are due to what we call the 'regular' (see section 8) parts of |QFS| 2 and |QFS| 2 , where the gluon is emitted outside the jet: dσ |QFS| 2 ; out; reg dijet dp + j1 d 2 p j1 dp + j2 d 2 p j2 (9.26) and dσ The Dirac traces appearing in the above two expressions are easily calculated using the identities (A.5) and (A.7), and read: RI: instantaneous gluon emission We find for the contribution due to the instantaneous emission of a gluon: Sudakov logarithms in our NLO cross section have the right form and can be absorbed into CSS. It turns out that this first step is already highly non-trivial, and we will limit ourselves in this work to the Sudakov double logarithms at large N c . To double leading logarithmic (DLL) accuracy, the TMD factorization formula for our process γA → dijet + X is: with some slight abuse of notation, since the product is really at the integrand level, and where S A is the perturbative Sudakov factor which reads at LO and in the DLL approximation: In the above formulas, b − b is the transverse coordinate conjugate to k ⊥ . We will show in this section that our NLO cross section, in the correlation limit, takes the form dσ NLO γA→dijet+X corr. lim. 4) hence proving agreement between the CGC and the TMD calculations, at least to DLL accuracy. The Sudakov double logarithm comes, at least at leading N c , from soft-collinear gluon radiation just outside the jet, i.e. from the contributions V out; phase |QFS| 2 (eqs. (8.20) and (8.21) and its q ↔q counterpart. However, as was remarked very recently in [120], taken at face value the double logarithm in this contribution comes with the wrong sign (see sec. 10.2). As we will demonstrate in subsection 10.3, this wrong sign is compensated for by the mismatch between naive and kinematically-improved low-x resummation.
Finally, we should remark that in the seminal paper ref. [49], eq. 10.4 was already inferred from an analysis of Higgs hadroproduction combined with kinematical arguments. While physically insightful, the approach to the dijet case in [49] has some limitations. The calculation presented in this section relies on the complete NLO dijet calculation instead, which requires a systematic treatment of the low-x resummation. As such, our treatment constitutes an important step towards reaching full NLO accuracy in the back-to-back limit.

Leading order
Introducing the vector sum of the transverse momenta of the outgoing (anti)quark: one can rewrite the transverse coordinates in function of b and r which are conjugate to k ⊥ and to P ⊥ , respectively (remember that z = p + 1 /q + andz = p + 2 /q + ): After this coordinate transform, the squared of the LO amplitude (2.23) becomes: × Tr Q 122 1 − s 12 − s 2 1 + 1 . (10.7) Taking the correlation limit P ⊥ k ⊥ then implies b, b r, r , which we will use to perform a Taylor expansion of the Wilson lines. At the lowest non-trivial order, the only nonzero contribution will come from the quadrupole since all O(r) corrections on dipoles evaluated in either two unprimed or two primed coordinates disappear: It is then easy to see that: up to higher orders in the Taylor expansion. eq. (10.7)) then simplifies to: (10.10) In the thus obtained Wilson-line structure, one can recognize the so-called 'hadron correlator' which is parameterized by the unpolarized F WW and linearly-polarized H WW Weizsäcker-Williams gluon TMD ( [61] (10.11) Finally, using eq. (2.1)), we end up with (note that at LO and at the cross section level, the parton momenta can be identified with jet momenta): which is the same expression as one would find in a leading-order TMD factorization approach [60]. Note that at the present lowest perturbative order, one needs an extra scale such as a nonzero photon virtuality Q 2 or heavy-quark mass to be sensitive to the linearly-polarized gluon TMD H WW [63,64,67].

Sudakov double logs in the NLO cross section
The contributions V out; phase |QFS| 2 and V out; phase |QFS| 2 (see eqs. (8.20), (8.21) and (8.25)) factorize from the LO cross section (at integrand level), and can produce double large logarithms. In the correlation limit, since p j1 → P ⊥ and and V out; phase has the same expression in the correlation limit. For simplicity let us restrict ourselves to the case in which the azimuthal angle of P ⊥ is integrated over. Then, we have in its integrand. These logarithms are indeed of the Sudakov type. However, the coefficient of the double log is positive, whereas the total Sudakov double log term should be negative. The only other class of diagrams in our calculation which can give a Sudakov double logarithm are the diagrams involving soft divergences and which give a contribution with the same Wilson-line operator as at LO. In the final-state exchange diagrams discussed in section 5, such as GEFS or the real interference between QFS and QFS, the leading-N c operators do not include quadrupoles (and drop in the correlation limit), and only subleading N c terms are proportional to the color operator appearing at LO.
In the case of high-energy logs in section 6, there is a cancelation of the subleading-N c terms between the C F terms from |QFS| 2 and |QFS| 2 , the final state diagrams discussed in section 5, and the real initial state emission diagrams |QSW + QSW| 2 . However, since our projectile, a photon, does not carry color charge, the initial state emission diagrams |QSW + QSW| 2 cannot have soft divergences. Hence, the cancelation of subleading N c terms occurring for high-energy logarithms cannot fully extend to the case of Sudakov double logarithms.
Due to the complexity of the diagrams from section 5, we are not attempting in this study to calculate their possible subleading-N c contribution to Sudakov logarithms, and we stay instead at large-N c . Hence, we can actually replace C F by N c /2 in eq. (10.15) at this accuracy.
In our calculation of the NLO dijet cross section in the correlation limit, we thus obtain Sudakov double logarithms with the expected coefficient at large N c but the wrong sign. We will now discuss the effect on this result from the kinematical improvement of the high-energy leading log resummation.

Sudakov double logs from the mismatch of naive and kinematically consistent low-x resummation
In order to perform the subtraction and resummation of low-x (or high-energy) logs for the NLO dijet cross section, we have used the simplest scheme for the JIMWLK resummation. As explained in section 6.1, in this scheme the evolution takes place along the p + axis, from a cutoff k + min to a factorization scale k + f . The JIMWLK evolution then resums multiple gluon emissions within this interval, which are strongly ordered in p + only.
However, such a simple scheme for the low-x resummation is known to fail beyond lowx leading logarithmic accuracy. Indeed, such a scheme amounts to an approximation of infinite collision energy √ s in order to simplify the kinematics, which then leads to serious issues like NLO cross sections which can go negative and unstable low-x evolution equations at next-to-leading log accuracy. The main ingredient for a resolution of these issues is an improvement of the kinematics in the high-energy evolution equation [86][87][88][89][90][91][92][93][94], also known as consistency constraint or kinematical constraint, which can be also understood as an all order resummation of leading collinear logarithms within the high-energy evolution equation, either BFKL [86][87][88][89][90] or BK [90][91][92][93]. The case of the JIMWLK equation has been much less studied so far, due to its complexity, although the same issue is present and can be addressed following a similar strategy [94].
The main reason for this failure is that the strong p + ordering of successive emissions is necessary but not sufficient for large high-energy leading logs to arise in higher-order calculations for a large but finite collision energy √ s. Instead, the necessary and sufficient condition is that successive emissions should be strongly ordered both in p + and in p − , in opposite directions. The JIMWLK equation in the scheme that we have used thus leads to an oversubtraction of high-energy logarithms, in the domain satisfying the p + ordering but violating the p − ordering. Note that the situation would have been similar if we had chosen a scheme for JIMWLK based on p − ordering only, except that the oversubtraction would have happened in a different kinematic domain.
The contribution to be subtracted from the NLO correction and resummed into the LL evolution of the LO term should thus be defined with two conditions, which require two factorization scales. Schematically, in addition to the condition p + 3 < k + f for the gluon of momentum p 3 to participate to the high-energy evolution, one should also include a condition of the type p − 3 > k − f . Note that both factorization scales k + f and k − f are used to specify the endpoint of the high-energy evolution of the target. They should thus be chosen in relation to the kinematics of the observable (the dijet in our case), and be independent of the kinematics of the target and of the collision energy √ s. First, k + f should be smaller or equal to the + momenta of both jets. A natural choice for k + f is then: The other factorization scale, k − f , should be a typical − momentum scale set by the dijet final state. Using the dijet mass squared (10.17) and q + , a natural choice for k − f can be written as: The BK and JIMWLK evolution equation are usually written in mixed space, with + momenta and transverse positions. Hence, one has no direct access to the p − of the gluon, which complicates the imposition of a condition of the type p − 3 > k − f . In practice, one is led to build a proxy for the p − of the gluon out of the available mixed-space variables, so that k − f becomes a lower bound on a combination of + momenta and transverse positions involved in the gluon emission [90][91][92][93][94]. In the present section, we do not need to enter into such technicalities. We can stay at a quite schematic level in order to understand the interplay between Sudakov double logs and the kinematical improvement of the JIMWLK equation.
Sudakov logarithms are known to exponentiate, and thus to appear at higher orders in terms containing the same Wilson-line operator (or the same parton distribution) as the LO term. Hence, we will now focus on the part in the JIMWLK evolution that involves this same LO operator. From the results of section 6, in our naive scheme for JIMWLK, at the cross section level this contribution amounts to multiplying the LO cross section with the factor In other sections, we have obtained the NLO dijet cross section using a naive scheme for high-energy log subtraction (and resummation). Using instead a kinematically improved scheme for this subtraction, we would obtain the same result for the NLO cross section, plus an extra term: the difference between the naive and the kinematically improved integrated JIMWLK evolutions of the LO cross section. The contribution to this difference which is proportional to the LO operator would then be with the term 1 in the bracket corresponding to the naive evolution with p + ordering only, whereas the theta function implements the condition (10.20), meaning imposing as well the p − ordering of the gluons participating to the evolution with respect to the k − f scale set by the dijet kinematics. Making all the bounds explicit as theta functions, the expression (10.21) becomes:

(10.22)
In terms of p + 3 , there is one upper but two lower bounds. k + min was first introduced as a cutoff to regulate the p + 3 -integral. In the high-energy resummation, k + min also plays the role of the starting point of the evolution of the target, and allows us to implement the dependence of the cross section on the energy of the collision as k + min ∝ 1/s. According to both of these interpretations for k + min , the lower bound p + 3 > k + min is less restrictive than the other one in the expression (10.22), and can be dropped. Then, eq.   for |P ⊥ ||b − b| 1. This is a Sudakov double log type term, this time with the expected negative sign. Combining this contribution together with the leading N c term from eq. (10.15), one finally obtains the following total double logarithmic contribution − α s N c 4π ln 2 P 2 ⊥ (b − b) 2 c 0 2 (10.25) at leading-N c , which is indeed the expected Sudakov double log term from eq. 10.4. Several remarks are in order: • The calculation outlined in this subsection provides an extra motivation for the kinematical improvement of high-energy evolution equations like JIMWLK: it allows one to obtain the correct Sudakov double logarithms as a leftover in the NLO cross section after high-energy resummation. The situation can be summarized as follows. In the naive scheme for JIMWLK resummation based on p + ordering only, any gluon radiation with smaller p + than the dijet scale k + f is treated as part of the evolution of the target. By contrast, including kinematical improvement allows one to split such small-p + gluon radiation into two parts: the true contribution to the evolution of the target with smaller p + but larger p − than the dijet, and the soft radiation with both p + and p − smaller than the dijet. Sudakov logarithms originate from the soft regime only, which is thus distinct from the true regime contributing to the evolution of the target. But the naive scheme in p + without kinematical improvement for the target evolution misses this fact, and leads to oversubtracting high-energy logs out of the soft regime. This motivates future studies in order to understand the practical implementation of the proposal of ref. [94] for the kinematical improvement of JIMWLK, or to construct other prescriptions.
• We have focused on the leading-N c contribution only. Subleading-N c terms cancel in the naive version of the BK and JIMWLK equation, but this cancelation can be broken when kinematical improvement is included, as can be seen from ref. [91] in the BK case. Hence, at this stage, we have no control on a possible subleading-N c correction to the coefficient of the Sudakov double log.
• If we had used a scheme for JIMWLK resummation based on p − ordering only, the situation would have been more favorable before kinematical improvement. In that case, only gluon radiation with larger p − than the dijet (and any p + ) would have been treated as part of the evolution of the target. This would not overlap with the soft regime characterized by both smaller p + and p − than the dijet. Hence, Sudakov logs should be obtained correctly in that case even without kinematical improvement of JIMWLK, since the oversubtraction of high-energy logs happens in the regime of larger p + and larger p − than the dijet, which is well separated from the soft regime. A formulation of the BK equation as evolution along p − was proposed in ref. [93]. However, it is crucially based on specific properties of BK, and for the moment such a scheme does not exist for JIMWLK.

Beyond the double leading logarithmic approximation
In principle, having obtained the full NLO cross section in general kinematics, one should be able to extend the notion of the correlation limit and write down a TMD-factorized NLO cross section for back-to-back dijets. All virtual diagrams would contribute since they preserve the leading-order kinematics. The most important real NLO corrections stem from gluon emissions inside or close to the jets, yielding Sudakov double-and single logarithms as well as finite terms. A first step is to extend the calculation above to the Sudakov single logarithms, which is left for future work. In the second step, where all the virtual diagrams need to be analyzed to obtain the finite contributions to the NLO cross section, we immediately encounter some difficulties. Indeed, TMD factorization is obvious for all diagrams with initial-or final-state loop corrections (i.e. IS + UV + FSUV, GEFS, IFS), at least in the sense that up to power corrections one can write the Wilson lines as the WW gluon TMDs (10.11) which decouple from the hard part. Note that the rationale for this power expansion in our approach comes from the phases e iP ⊥ ·r e ik ⊥ ·b which imply that r b when P ⊥ k ⊥ . However, in all virtual graphs where the gluon scatters off the shockwave (SESW, GESW, ISW) this procedure is compromised due to the phase e iP ⊥ ·(r+ξx 13 ) e ik ⊥ ·b , which now enforces |r + ξx 13 27) and would bring the virtual NLO contributions in a TMD-factorized form. Unfortunately, it is not clear to us whether such an expansion can be justified.

Conclusions
Making use of the dipole picture of the CGC effective theory and of LCPT, we have calculated the cross section for the inclusive production of two jets in the scattering of a real photon with a target proton or nucleus at low x. The computation was performed at next-to-leading order in α s , while resumming the multiple rescatterings of the partons off the semiclassical gluon fields in the target to all orders in the eikonal approximation.
Using dimensional regularization, we explicitly showed the cancellation of ultraviolet singularities between different virtual NLO corrections on the amplitude level. Likewise, we demonstrated how both soft and collinear divergences, which appear in final-state radiation, cancel through an intricate interplay between the jet definition and certain virtual diagrams. Moreover, we regularized rapidity divergences with the standard cutoff method and absorbed the resulting large logarithms into the JIMWLK equations. This hybrid scheme, i.e. the dimensional regularization of UV and soft-collinear divergences combined with a cutoff for rapidity divergences, is commonly used in higher-order CGC calculations. However, it has the drawback that it obfuscates the distinction between 'genuine' soft divergences on the one hand, and rapidity divergences on the other. The former are typically associated with initial-and final state radiation, while the latter are the hallmark of low-x physics and related to the highly boosted target. The fact that, as we show in section 8, not only all soft and collinear singularities cancel in the final-state, but the only large logarithms left are those removed by JIMWLK, is a powerful confirmation of the consistency of this scheme. After having obtained the NLO dijet cross section, we have explored the back-to-back limit to investigate whether our result could be cast in a form consistent with TMD factorization. At leading order, the overlap of the CGC and the TMD frameworks for processes such as this has been demonstrated already some time ago [59,60]. However, a full nextto-leading order matching is much more involved, since it constitutes an analysis of all Sudakov double and single logarithms as well as finite NLO contributions. A first demonstration that Sudakov double logarithms arise in the hadroproduction of a Higgs boson at low-x was performed in [49], where also the precise form of these double logs in dijet production was inferred based on kinematical arguments. In this work, we revisited the analysis of the Sudakov double logarithms in the dijet case, this time based on the full NLO calculation. We argue that a kinematical improvement of the JIMWLK resummation is crucial to obtain the correct result for the Sudakov double logarithms. We have also identified a class of virtual contributions that, at least at first sight, break TMD factorization on the finite or potentially single-logarithmic level. Before drawing definite conclusions a more thorough study of the correlation limit of our process is needed, which is left for future work. The inclusive dijet electroproduction process has been studied within TMD factorization in refs. [121,122].
While this work was in progress, the NLO calculation of the inclusive dijet electroproduction cross section appeared in ref. [83]. This calculation was performed in a covariant formulation of the CGC rather than in LCPT, and using a different UV subtraction scheme. In the photoproduction limit Q 2 → 0, our cross section and the γ * T A → dijetX cross section in [83] should coincide. The largest difference between both results stems from the treatment of the jet. Indeed, the final result eq. 7.16 of [83] for the dijet cross section is still sensitive to both single-and double logarithms in the rapidity renormalization scale z f = k + f /q + (or in the rapidity cutoff k + min if the JIMWLK subtraction is performed after the application of the jet algorithm.) Since these logarithms cannot be absorbed into JIMWLK and have a soft origin, not a rapidity one, they are unphysical and still need to cancel with soft gluon emission outside the jet, as demonstrated in section 8. Other differences between our results are relatively minor and mainly related to the precise + momentum due to the Dirac traces. In appendix A.2, we cast our partonic cross section in the same notations and conventions as [83] to facilitate a detailed comparison.