NLO impact factor for inclusive photon + dijet production in e + A DIS at small x

We compute the next-to-leading order (NLO) impact factor for inclusive photon +dijet production in electron-nucleus (e+A) deeply inelastic scattering (DIS) at small x. An important ingredient in our computation is the simple structure of “shock wave” fermion and gluon propagators. This allows one to employ standard momentum space Feynman diagram techniques for higher order computations in the Regge limit of fixed Q ΛQCD and x → 0. Our computations in the Color Glass Condensate (CGC) effective field theory include the resummation of all-twist power corrections Qs/Q , where Qs is the saturation scale in the nucleus. We discuss the structure of ultraviolet, collinear and soft divergences in the CGC, and extract the leading logs in x; the structure of the corresponding rapidity divergences gives a nontrivial first principles derivation of the JIMWLK renormalization group evolution equation for multiparton lightlike Wilson line correlators. Explicit expressions are given for the x-independent O(αs) contributions that constitute the NLO impact factor. These results, combined with extant results on NLO JIMWLK evolution, provide the ingredients to compute the inclusive photon + dijet cross-section at small x to O(α s ln(x)). First results for the NLO impact factor in inclusive dijet production are recovered in the soft photon limit. A byproduct of our computation is the LO photon+ 3 jet (quark-antiquark-gluon) cross-section.

An important discovery of the electron-proton (e + p) deep inelastic scattering (DIS) experiments at HERA was the rapid growth of the gluon distribution with decreasing Bjorken x, for fixed large momentum transfer squared Q 2 . This demonstrated that the proton wavefunction in the corresponding high energy Regge limit is dominated by Fock state configurations containing large numbers of gluons. Their number grows via bremsstrahlung with increasing energy or decreasing x. It was conjectured that in this Regge limit repulsive gluon recombination and screening effects [1,2] conspire to slow down the growth of cross-sections. A remarkable consequence is that a semi-hard saturation scale Q s (x) is generated dynamically by these competing many-body effects.
If the saturation scale is large compared to the intrinsic QCD scale, asymptotic freedom suggests that the coupling α s (Q s ) 1; this allows weak coupling effective field theory techniques to be employed in systematically computing cross-sections in a regime of QCD where field strengths are nonperturbatively large. The physics of this nonlinear regime of QCD can be quantified in a classical effective field theory (EFT) called the Color Glass Condensate (CGC) [3][4][5][6][7][8][9] which implements a Born-Oppenheimer separation of the relevant degrees of freedom into static color sources at large x and gauge fields at small x.
A further important element in the EFT is that the separation in x between color sources and fields satisfies a renormalization group (RG) equation as the scale separation is evolved towards smaller x. This can be understood in a Wilsonian picture wherein, with each small change δx in the scale separation, the dynamical gauge degrees of freedom within δx are absorbed into the static light cone color sources in the classical EFT at the lower x − δx scale. The RG equation correspondingly describes the change in a nonperturbative weight functional describing the distribution of color sources, from large x towards smaller x, and efficiently resums simultaneously logarithms in α s ln(1/x) and power corrections Q 2 s (x)/Q 2 when they become large. To leading logarithmic accuracy in x, this RG equation is the JIMWLK equation [10][11][12][13] with a corresponding JIMWLK Hamiltonian [14] that contains all the relevant information regarding the rapidity evolution of n-point Wilson line correlators.
A number of well-known results are obtained as limits of the JIMWLK Hamiltonian. In the limit of large number of colors N c , and for large atomic nuclei with A 1, the JIMWLK equation for the simplest "dipole" correlator of light-like Wilson lines describing fully inclusive DIS is the Balitsky-Kovchegov (BK) equation [15,16]. In the leading twist limit, where Q 2 s (x)/Q 2 1, this reduces to the BFKL equation [17,18]. The latter was first derived by explicit computation of perturbative QCD Feynman diagrams in Regge asymptotics.
The CGC EFT has been applied to compute a large number of final states in both DIS and in hadron-hadron collisions. An excellent introductory review of the formalism and applications can be found in [9]. An important test of the framework will be its ability to make predictions of high accuracy that can be compared to experiment. The ideal experimental conditions for such tests require access to large Q 2 and small x. Further, since the saturation scale is enhanced in nuclei as Q 2 s,A ∼ A 1/3 Q 2 s,proton [1,3,8] , heavy nuclear targets are preferred. In principle, these conditions are achieved in proton-nucleus (p + A) collisions at the LHC. However large final state interactions may be present in these experiments that will complicate interpretations of the data [19]. These considerations provide a major motivation for future Electron-Ion Collider (EIC) experiments [20][21][22].
Such experiments at an EIC are envisaged to have much higher luminosities than were available at HERA; this, in combination with the nuclear beams, greatly enhances the possibility that DIS e + A collider measurements may uncover definitive evidence for gluon saturation. These precision DIS experiments demand higher order computations in Regge asymptotics in close analogy to how higher order computations in the Bjorken limit provided powerful tests of perturbative QCD. To further the analogy, just as one computes universal splitting functions and process dependent coefficient functions in the Bjorken limit of DIS, one needs to compute both universal multiparton lightlike correlators (or equivalently, as we shall discuss, their generating weight functional describing color charge distributions) and process dependent impact factors to higher order accuracy in the Regge limit.
The next-to-leading order (NLO) computation of the BFKL kernel has now been available for over twenty years [23,24]. Subtleties in the treatment of singularities in the NLO BFKL kernel were noticed shortly after [25,26]-for a recent comprehensive discussion, see [27]. Computations of higher order corrections to the BK equations followed in short order [28][29][30] and specific prescriptions for the running coupling following from these studies were implemented [31] in phenomenological studies. More recently, NLO computations of multi-point Wilson line correlators, or equivalently, the NLO JIMWLK Hamiltonian have become available [32][33][34][35][36] and next-to-next-to-leading order (NNLO) computations of BK and JIMWLK are underway [37,38].
Because of the absence of final state interactions, isolated photons are clean probes of this strongly correlated gluonic matter. Several computations [39][40][41][42][43] of inclusive photon production have been performed in the CGC framework in the context of proton-nucleus (p + A) collisions. In a recent paper [44], henceforth referred to as Paper I, we reported on a first CGC computation of the leading order differential cross-section for inclusive prompt photon production in conjunction with two jets in electron-nucleus (e + A) DIS at small x. This process has clean initial and final states and is the simplest non-trivial process besides fully inclusive DIS to study the physics of gluon saturation in e + A collisions. This computation provides more differential phase space distributions, thereby going beyond existing small x computations [45][46][47][48][49][50][51][52][53] on the total cross-section for fully inclusive DIS; exceptions are the NLO differential cross-section computations by Boussarie et al. [54][55][56], albeit for diffractive DIS.
In Paper I, we showed explicitly how in the limit of the final state photon momentum k γ → 0, the Low-Burnett-Kroll soft photon theorem [57][58][59] allows one to recover existing results for inclusive dijet production in DIS [60]. In the leading twist limit, we also obtained the k ⊥ and collinear factorized expressions which match the dominant NLO small x perturbative QCD (pQCD) contributions. In particular, our result in the collinear limit is directly proportional to the nuclear gluon distribution at small x.
For precision physics in pQCD, it is essential to go beyond leading order descriptions for quantitative studies of data. NLO computations will be especially important for the discovery and characterization of gluon saturation in e+A DIS where its effects are anticipated to be larger than in e + p DIS 1 . As we observed in our LO photon+dijet computation, there are novel quadrupole gauge invariant correlators of lightlike Wilson lines that appear, whose energy evolution, in addition to the dipole correlators measured in fully inclusive DIS, will be a sensitive test of JIMWLK evolution. We will show how such correlators, in combination with the dipole correlators, violate the soft gluon theorem. A quantitative understanding of this violation can provide deeper insight into the infrared structure of QCD in the Regge limit [61].
Byproducts of our computation of the differential photon+dijet cross-section are the first NLO results for inclusive dijet, inclusive photon and photon+jet measurements at an EIC. Further, the NLO graphs for real gluon emission provide the complete LO results for the photon+3-jet (γ + qqg) and 3-jet qqg final states [62]. We will point to the steps necessary for extracting "numbers" from our computation; though much more computationally challenging than comparable computations in the highly developed collinear factorization pQCD framework, such a program is feasible and essential for precision physics at an EIC.
As we will discuss further shortly, all computations in the CGC EFT rely on a separation of scales between static color sources and dynamical gauge fields. Thus perturbative computations at small x in this framework are performed in a background of such static color sources and physical quantities are obtained by subsequent gauge invariant averaging over these sources. The first principles formalism in quantum field theory underlying such computations in strong background fields has been discussed previously [63][64][65][66][67]; in particular, [63] and [67] provide pedagogical discussions in complementary approaches.
We begin our discussion of the NLO DIS photon+dijet computation with the starting point of all CGC computations, the classical Yang-Mills equations, Here the covariant derivative D µ = ∂ µ − igA µ , g represents the QCD gauge coupling, and ρ A (x − , x ⊥ ) ≈ ρ A (x ⊥ )δ(x − ) represents the color charge density of large x static sources for the small x dynamical fields A µ . The delta-function in the color charge density denotes that we are working in a frame where the nucleus is moving in the positive z-direction at nearly the speed of light with large light cone longitudinal momentum P + N → ∞. (See Appendix A for details of the conventions adopted in this work.) In addition we will choose the frame in which the virtual photon has a large longitudinal momentum q − and transverse momentum q ⊥ = 0.
The solution of the classical equations in Lorenz gauge ∂ µ A µ = 0 is given by where Λ is an infrared cutoff that is necessary to invert the Laplace equation −∇ 2 ⊥ A + cl = gρ A in two dimensions. This solution to the Yang-Mills equations in Lorenz gauge is simply related to the solution in the light cone gaugeÃ + = 0, where one obtains likewise thatÃ − cl = 0 andÃ i cl = i g U ∂ i U † , and denotes the adjoint Wilson line expressed in terms of the the large x static color source densities via Eq. (2). Note that T a , a = 1, · · · , 8 are the generators of color SU (3) in the adjoint representation. These Wilson lines, and their counterpartsŨ in the fundamental representation represent respectively, the path ordered phase acquired by a gluon and a quark in their eikonal multiple scattering off the classical background field of a nucleus. The Wilson lineŨ is obtained by replacing the adjoint generators in Eq. (3) with the fundamental generators: T a → t a . In the LO photon+dijet cross-section of Paper I, the virtual photon fluctuates into a quark-antiquark pair that multiple scatters off the classical background field of the nucleus. In the Feynman diagram computations of this LO process, the Wilson lines are incorporated in the momentum space structure of the dressed quark propagators. At NLO, there are real and virtual gluon contributions to the leading order process. The corresponding gluon propagators are also dressed by multiple scattering off the classical background field of the nucleus.
The structure of the dressed quark and gluon propagators in the classical background field, is particularly simple in the "wrong" light cone (LC) gauge A − = 0. As indicated by Eq. (2), this gauge shares the same classical field solution with Lorenz gauge. However, unlike the LC gauge A + = 0 for P + N → ∞, it does not provide a simple physical interpretation of parton distribution functions. Any concern in this regard is however far outweighed by the advantages provided by the simple forms of the dressed propagators that were computed previously in Refs. [5,15,[68][69][70][71]. The effective vertices for these dressed quark and gluon "shockwave" propagators, incorporating the fundamental and adjoint Wilson lines respectively, are shown in Fig. 1.
As discussed in Paper I, these effective vertices are identical to the quark-quark-reggeon and gluon-gluon-reggeon effective vertices [36,[72][73][74][75][76] in Lipatov's reggeon field theory [77]. Another salient feature in the expressions below is that we do not subtract the unit matrix in the expansion of Wilson lines. This allows for the possibility that the quarks and gluons do not scatter in addition to all their possible multiple scatterings encoded in the higher terms in the Wilson line expansion. Consequently, we draw Feynman diagrams with all dressed fermion and gluon propagators and for each such kinematically allowed process, we only need to subtract the "no scattering" contribution (obtained by puttingŨ and U 's to unity) to get the physical amplitude. As we will see, this significantly aids in the NLO computation where the number of contributing processes is large.
Note further that since our computations do not employ light-front perturbation theory like many of the NLO computations in the literature, and are carried out entirely in momentum space (also unlike many computations), they also provide a useful cross-check on extant NLO results on fully inclusive DIS. The techniques employed here may also provide a pathway to carrying out higher order computations to NNLO in the Regge limit.
To proceed with our NLO computation, it is important to elaborate further on the RG procedure for resummation of large logarithms in x, specifically with regard to how it applies to the inclusive photon+dijet computation of interest. As noted, the Wilsonian RG ideology on which the CGC EFT is based naturally involves a cutoff scale in rapidity or longitudinal momentum separating the soft and hard partons in a hadron/nucleus. At LO, this scale Λ + 0 (or rapidity Y 0 = ln(Λ + beam /Λ + 0 )) is arbitrary and the fast or valence modes with longitudinal momenta k + Λ + 0 are represented by the stochastic color charge density, ρ A (x − , x ⊥ ). A gauge invariant weight functional W Λ + 0 (Y0) [ρ A ] describes the probability density corresponding to this charge density. As also noted, the soft modes are represented by classical color fields that are solutions of the classical Yang-Mills equations with appropriate gauge fixing conditions. As we boost the nucleus towards the small x scale of interest (or towards higher energies), modes that were below the cutoff Λ + 0 start contributing to the scattering process. We therefore have to consider quantum effects induced by these "semi-fast" gluons [12] (see Fig. 2) which can be defined as the nearly on-shell fluctuations with momenta deeply inside the strip Λ + 1 (= bΛ + 0 ) |l + | Λ + 0 or conversely, energies in the range Here Q 2 0 and x 0 are respectively the virtuality of the nucleus and the Bjorken-x at the initial scale. Further, Q 2 = −2q + q − is the fixed virtuality of the exchanged virtual photon in DIS (in Regge kinematics) and x is the Bjorken-x of interest determined by the kinematics of the process. The effect of integrating out these fluctuations manifests itself in the appearance of large logarithms ln(1/b) which for α S ln(1/b) ∼ 1 must be resummed to all orders in α S . Schematic illustration of sources and fields in the CGC effective theory in terms of the cutoff scale Λ + 0 , or equivalently the rapidity Y0. At NLO, contributions from field modes in the range Λ + 1 < k + < Λ + 0 , such that αS ln(Λ + 0 /Λ + 1 ) ≡ αS∆Y < 1, are integrated out and absorbed into the source densities at the scale Λ + 1 . This self-similar renormalization group (RG) pattern is repeated successively generating the JIMWLK RG equation for the source densities.
Denoting the differential cross-section for inclusive photon plus dijet production by dσ for simplicity, we can write down its expectation value in the CGC EFT at LO as [44] The r.h.s represents the fact that the LO cross-section is first computed for a fixed distribution of color sources ρ A with Λ − < Λ − 0 . This object dσ LO [ρ A ] shown in Fig. 3, is computed using standard techniques in perturbative QCD albeit, as noted, with the modified propagators listed in Fig. 1, wherein the dependence on ρ A enters (at LO) via the fundamental Wilson lineŨ . 3. Processes contributing to the leading order amplitude and hence the cross-section dσLO in Eq. 5. The other two diagrams are obtained simply by interchanging the quark and antiquark lines.
The process independent weight functional W Λ − 0 [ρ A ] is a nonperturbative object that contains fundamental information about n-body correlations amongst gluons at the initial scale Λ − 0 . It can be understood as representing large x diagonal elements of the density matrix of QCD in the Regge limit; a recent discussion of W [ρ A ], and generalizations thereof, can be found in [78]. At NLO (≡ O(α S )) in the CGC power counting, we have to account for quantum fluctuations of both the quarkantiquark dipole as well as the wavefunction of the nuclear target. In the EFT language, these are distinguished by the magnitude of the LC momentum of the gluon modes relative to the initial scale Λ − 0 . As shown in Fig. 4, the modes with l − < Λ − 0 (which we shall denote as NLO:1) can be interpreted as contributions from Fock states dressing the target wavefunction. The contribution to the cross-section for processes of this kind can in general be written as where dσ NLO:1 [ρ A ], for a fixed configuration of ρ A , is comprised of nontrivial combinations of Dirac traces and Fourier transforms of color traces over products of the Wilson linesŨ and U . From these contributions, we are interested in collecting only pieces that contain large logarithms in Λ − 0 . These can be written as where Λ − 1 is the scale to which the target gluon modes are evolved. Here H LO represents the JIMWLK Hamiltonian [10][11][12][13][14]; its explicit form will be discussed later in the paper. For our purposes here, it suffices to note that H LO dσ LO is of order α S dσ LO [ρ A ]. Combining the above contribution with the LO result in Eq. (5), and using the Hermiticity of W with respect to the functional integration over ρ A , we can write the result as Further redefining and thereby absorbing the effects of the semi-fast gluons in terms of a modification of the probability distribution of the fast color sources, one obtains the leading log in x JIMWLK equation To derive this result, we employed the essential RG philosophy that the observable on the l.h.s of Eq. 8 must be independent of the arbitrary scale Λ − 0 separating the static color sources ρ A from the dynamical gauge fields. Replacing the expression in curly brackets in Eq. (8) by the r.h.s of Eq. (9) is equivalent to summing the leading logarithmic terms α S ln(1/x) to all orders in perturbation theory. We will henceforth label the weight functional that satisfies Eq. (10) as W LLx [ρ A ].
One also has NLO contributions from gluon modes with l − > Λ − 0 corresponding to quantum fluctuations of the dipole projectile. As shown in Fig. 5, these include real gluon emission and virtual gluon exchange processes; the filled blobs represent the dressed gluon propagators allowing for the possibility that the gluon can scatter off the background classical field of the nucleus. We will refer to contributions from the quantum fluctuations of the projectile as NLO:2 contributions to distinguish these from the NLO:1 quantum fluctuations of the target below the scale Λ − 0 .
Representative NLO contributions from gluon modes in the projectile with l − > Λ − 0 . Both real and virtual emission diagrams are shown. For the latter case, we have to consider interference diagrams with LO processes described in [44].
As in any loop computation, the intermediate steps of our calculation will contain soft, collinear and ultraviolet (UV) singularities depending on the region of phase space of the gluon that we are integrating over. In the virtual graphs, UV divergences appear from integrals over the transverse momentum of the gluon in the loop; these are isolated using dimensional regularization in d = 2 − dimensions. At this NLO order of quantum fluctuations of the projectile quark-antiquark pair, all UV divergences must vanish or cancel without the necessity of renormalizing the parameters of the EFT. This is because we will work the limit of massless quarks and there is no running of the QCD coupling constant in the projectile wavefunction at this order in the CGC power counting.
The small x divergences arise from integrating over the quantum fluctuations induced by "slow" or "semi-fast" gluons with '-' longitudinal momenta that are small relative to the large q − momentum of the virtual photon. These divergences are regulated by imposing a cutoff at the initial scale, Λ − 0 of the evolution which is defined in Eq. 4. The resulting logarithms in Λ − 0 (or equivalently x) are absorbed into small x renormalization group evolution of the weight functional W [ρ A ] as shown in Eq. 10. At higher orders, it may be necessary to employ more sophisticated regularization schemes for these rapidity divergences as well [79].
For gluon emission diagrams, in addition to small x divergences, there are also singularities that arise from the region of phase space where the unscattered gluon is soft or collinear to the quark or antiquark. In particular, there are residual collinear divergences that survive after real and virtual contributions are combined. These divergences are absorbed into the evolution of fragmentation functions. Conversely, we can regulate the phase space integration over final states by promoting partons to jets where the latter are defined using a cone algorithm [80,81]. We will show explicitly in the limit of small jet cone size [82][83][84][85] that collinear divergences between real and virtual graphs cancel completely enabling the extraction of the dominant contributions towards the jet cross-section.
With all the divergences in the (NLO : 2) quantum fluctuations of the virtual photon projectile accounted for, one can write the infrared (IR) safe jet cross-section as These NLO contributions (shown in Fig. 5) can be broken up into two pieces. The first piece is obtained by taking the "slow" gluon limit, l − → 0 and is identical (δσ jet NLO:2 = δσ jet NLO:1 ) to the expression in Eq. 7 at the momentum scale Λ − 0 . We will show this explicitly later in the paper. Specifically, this matching corresponds to a first principles derivation of leading log JIMWLK evolution for a nontrivial final state of the projectile. While the NLO:2 derivation represents the slow IR limit of the projectile, the RG corresponds to matching it to quantum fluctuations at their fast scale Λ − 0 in the target.

8
The second term on the r.h.s of Eq. 11 (dσ jet;finite NLO [ρ A ]) contains genuine α S suppressed (without logs in x) contributions to the differential cross-section from real and virtual graphs. For the latter, it is possible to deduce analytical expressions because the divergent structures can be isolated at the level of the amplitude. In contrast, divergent structures in the real NLO contributions are manifest only at the level of the squared amplitude and obtaining analytical expressions for the finite terms is challenging. However they can be evaluated numerically using the fact that rapidity divergences can be isolated in the slow gluon limit; these can then be subtracted from the squared amplitudes (using a numerical cutoff procedure) to obtain the desired finite pieces.
Further, by replacing parton momenta in these contributions with those of jets (using a jet algorithm) gets rid of the remaining collinear divergences. The finite contributions dσ jet;finite NLO [ρ A ] that we will compute explicitly in this paper are, in the language of Regge theory, the NLO "impact factor" corrections to the LO impact factor dσ jet LO [ρ A ]. Computing this process-dependent NLO impact factor dσ jet;finite NLO [ρ A ] for photon+dijet production is important because it allows us to go one step further in precision and consider relevant (two loop) NNLO contributions to the cross-section that have terms proportional to α 2 S ln(Λ − 1 /Λ − 0 ). These contributions are effectively of NLO magnitude if α S ln(Λ − 1 /Λ − 0 ) ∼ 1. Diagrams corresponding to a two loop fluctuation of the target are shown in Fig. 6. In the class of such two loop diagrams, there are contributions of order which are included in the leading log JIMWLK resummation, as represented by W LLx [ρ A ]. There are also contributions from two-loop QCD diagrams proportional to α 2 S alone (without leading logs in x) but these are suppressed at the desired accuracy of our problem. We will consider here only those two loop contributions in Fig. 6 that contain next-to-leading logarithms in x (NLLx) contributions to the result in Eq. 8. This in turn gives us the LO+NLLx result which can be expressed in terms of a modified weight-functional as where and the NLO JIMWLK Hamiltonian H NLO [32][33][34][35][36] is of order α 2 S . 6. NNLO diagrams that contain αS magnitude contributions for αSln(Λ − 1 /Λ − 0 ) ∼ 1. These diagrams contribute to the NLO JIMWLK kernel.
Representative NNLO diagrams whose leading logarithmic pieces are combined with the αS suppressed contributions to show the LLx evolution of the weight functional for the large-x color sources at NLO.
There are however a second class of two loop NNLO contributions of the sort shown in Fig. 7 which contain contributions that are parametrically of order α 2 S ln(Λ − 1 /Λ − 0 ). These correspond to one loop fluctuations of both the projectile and the target. From these processes, we have to extract the LLx pieces from the gluon fluctuations below the cut Λ − 0 and match them with the O(α S ) NLO "impact factor" expression in Eq. 11 to obtain As a result of our power counting in powers of α S and α S ln(Λ − 1 /Λ − 0 ), one can finally write the complete NLO result for the differential cross-section at NLLx accuracy by combining the results in Eqs. 12 and 14 respectively as As the in the above equation indicates, we can go one step further by promoting the second term on the r.h.s of the first equality from W LLx → W N LLx . This extends the scope of computation to O(α 3 S ln(1/x)) with the understanding that we will miss terms at that order of accuracy.
The weight functional W N LLx [ρ A ] in Eq. (13) can be obtained by adapting extant results for LO and NLO(≡O(α 2 S )) JIMWLK [33,35,36,86,87] and BK [28,29,32] evolution into our approach. Thus to obtain results for photon+dijet production up to O(α 3 S ln(Λ − 1 /Λ − 0 )) accuracy, it is sufficient to compute the NLO impact factor dσ jet;finite To go beyond this level of accuracy, we would have to compute fluctuations in the projectile involving the emission of two gluons, two loop virtual gluon processes, as well as interference diagrams of the genuine NLO processes shown in Fig. 5. This will be reserved for a future project.
The rest of the paper is organized as follows. In Sec. II, we shall briefly review the essential elements of the leading order (LO) computation of Paper I [44] and revisit some of the key results obtained there. In Sec. III, we will outline the structure of the NLO computation, structured for convenience in two subsections. In Sec. III A we categorize contributions to the NLO amplitudes from real gluon emission and virtual gluon exchange diagrams in terms of their color structure. We also provide an interactive flowchart (see Fig. 17) in this subsection furnished with hyperlinks that direct the reader to the final expressions later in the text for the various components constituting the real and virtual contributions to the NLO amplitude. In Sec. III B, we organize these contributions at the level of the squared amplitude in Table II in terms of common Wilson line structures. Doing so enables one to see cancellations of divergences in a transparent manner; this in turn facilitates the computation of the finite NLO impact factor in the photon+dijet inclusive cross-section.
Section IV contains a detailed computation of the amplitude for real gluon emission processes. These are discussed separately for the case when the real gluon either crosses or does not cross the nuclear "shock wave" using a representative diagram from each category. The final expression for the amplitude is given by Eq. 75. In section V, we describe in detail the computation of the amplitude for virtual gluon exchange processes; these are broadly classified into the topologies of self-energy and vertex corrections. Specifically, sections V A and V B deal with the amplitudes for self-energy graphs with dressed and free gluon propagators respectively. There are ultraviolet and rapidity singularities associated with these processes which are carefully isolated from the finite parts. For each class of diagrams, we use a representative graph to show the explicit computation. The results for the amplitudes from these processes are given respectively by Eqs. 79, 106 and 136. A similar exercise is performed in sections V C and V D respectively for the vertex correction processes with dressed and free gluon propagators. The final expressions for the amplitudes are given by Eqs. 139 and 164 for the case of dressed gluon propagators and by Eqs. 167 and 190 for the case in which the gluon does not cross the shock wave.
Section VI combines the results obtained in the earlier sections to obtain the final result for the principal goal of our study, the NLO impact factor for photon + dijet production in e + A DIS. We demonstrate here the cancellation of collinear divergences between real and virtual processes resulting in an infrared safe differential cross-section. To facilitate this, we introduce jet definitions and work in the approximation of a jet with small cone radius [82] to explicitly extract the collinearly divergent contributions from the squared amplitudes of real gluon emission graphs which contain the possibility of a gluon being collinear to the (anti) quark. We note that there is no Sudakov suppression of the cross-section because we have not imposed any kinematic constraints. Interestingly, we observe that (unlike the case of diffractive DIS [55]) the NLO cross-section does not factorize into the LO result and kinematic factors in the soft gluon limit. The implications of this result will be addressed in future work.
In Section VII, we take the slow gluon limit of our general expressions for the cross-sections and show that these provide a first principles derivation of the JIMWLK RG equation. While there exist several derivations of the JIMWLK equation in the literature going back to the original papers, many of these begin at the outset in the slow gluon limit. It is therefore interesting to see how the JIMWLK equation arises in the explicit computation of the nontrivial photon+dijet cross-section. This exercise also helps lay the groundwork for an independent derivation of the NLLx JIMWLK equation.
We will end this paper with a brief summary and outlook. With regard to the latter, an important next step is to provide quantitative predictions for measurements at a future EIC. These are significantly more challenging even though our use of the wrong light cone (A − = 0) gauge allows us to present our computations in a manner analogous to comparable NLO computations in collinear factorization computations. This is firstly because going away from the collinear limit introduces additional nontrivial integrals in the computations. Further, a quantitative computation of the dipole and quadrupole correlations is much more complex than their parton distribution (pdf) counterparts. This is unsurprising because the former contain a tremendous amount of information on the physics of many-body correlations in QCD that are not contained in the pdfs. Nevertheless, the technology to achieve the desired goal has advanced considerably to bring it within reach.
The principal results and conclusions of this paper are spelled out in an accompanying letter [88]. Appendices A through J supplement the material in the body of the paper. The notations and conventions used throughout the paper are summarized in Appendix A. In the computation of the amplitude for the various processes, we will encounter tensor integrals over transverse components of the gluon loop momenta. General expressions for these constituent integrals along with details for special cases are provided in Appendix B for both processes with gluon emissions and gluon loops. Appendix C contains detailed expressions for the amplitudes for gluon emission processes that are too cumbersome to include in the main text. Likewise, in Appendix D, we provide a detailed computation of the quark self-energy, which provides the template to compute the amplitudes of self-energy graphs where the gluon propagator is not dressed. The expression obtained in Eq. D1 for the gluon loop contribution is very general and can be straightforwardly used in any pQCD computation performed using light cone coordinates and in the light cone gauge. A similar computation for the virtual gluon corrections to the γqq and γ * qq vertices is provided in Appendices E and F.
Appendix G contains the rapidity divergent pieces, discussed in Sec. V C, for the vertex corrections with the dressed gluon propagator that are not provided in the main text. Similar expressions for the amplitudes with final state interactions (discussed in Sec. V D) are provided in Appendix H. The expressions for the finite pieces of the amplitudes are distributed over seven subsections in Appendix I. Finally, Appendix J provides a short proof of the sub-dominance of non-collinearly divergent contributions to the cross-section for real gluon emissions when we work in the limit of small jet cone radius.

II. GENERAL DEFINITIONS AND BRIEF REVIEW OF LO COMPUTATION
We will work in the light cone (LC) gauge A − = 0 throughout this computation. The highly energetic nucleus is considered to be right moving so that it has a large '+' component of LC momentum P + N . The virtual photon exchanged between the electron and nucleus is considered to be left moving and consequently has a large '−' component of LC momentum q − . The mass of the electron is neglected throughout the calculation.
Following the LO computation in [44], we can write the amplitude for inclusive photon+dijet production in DIS as where is obtained by index contraction with the propagator for the exchanged photon with momentum 2 q = (−Q 2 /2q − , q − , 0 ⊥ ). The amplitude for the hadronic subprocess is given by and is the quantity of interest. Here (k γ , λ) is the polarization vector for the outgoing photon. The 4-momentum 3 assignments are given in Table I and boldface letters denote 3-momentum vectors. We define the following ratios of the outgoing momenta to the dominant component q − of the incoming virtual photon momentum We will work in the limit of light quarks and neglect their masses in the present computation. Since we are dealing with prompt photon production in DIS at small x, the dominant contribution is indeed expected to involve light quarks. Squaring the expression for the amplitude in Eq. (16), and performing the necessary averaging and sum over electron spins and photon polarizations 4 , we can write The lepton tensor given by is identical to the one obtained in fully inclusive DIS. In the following, we will concentrate on obtaining the NLO contributions to the hadron tensor which is defined as The . . . in the above equation refers to the CGC averaging over all possible source charge configurations ρ A . From a first principles Quantum Field Theory perspective, this corresponds to the systematic computation of Feynman diagrams in the presence of static sources, and subsequently performing averages over the source distribution, as spelled out in [63][64][65], and references therein. For a generic operatorÔ this is quantified as [12,13] In this equation,Ô[ρ A ] is the quantum expectation value of the operator for a given charge configuration ρ A . One then performs the classical-statistical average ofÔ over all possible color charge configurations with the gauge invariant weight functional W Y [ρ A ] representing the distribution of the color charge configurations at a rapidity Y = ln(x/x 0 ) in the target. This double average is justified because the color charges ρ A are long-lived (or static) on the time scales corresponding to the (quantum) dynamics of the gauge fields. The functional dependence on ρ A enters the amplitude M µα through Wilson lines which are also the phase rotations in color space obtained by the quark and antiquark during their eikonal propagation along the light cone. We will see this more clearly when we present the structure of the amplitudes at LO and NLO in the upcoming discussion.
Since we wish our presentation to be self-contained, we will sketch here the LO contributions to the amplitude derived in [44]. At LO in the CGC power counting, there are four contributions to the amplitude; two of these 4 We use here the identity λ β (kγ , λ) * α (kγ , λ) = −g αβ + k α γ n β + k β γ n α k − γ , as the sum over outgoing photon polarizations. By virtue of the Ward identity, we can easily show that only terms proportional to g αβ contribute, thereby leading to Eq. 22. are shown below in Fig. 8 and the other two obtained by interchanging the quark and antiquark lines. As noted previously, an important ingredient in the computation is the simple form of the dressed quark propagator in the classical background field of the target nucleus [5,15,68,69,71]. In A − = 0 gauge, this can be expressed as [70] S ij (p, q) = S 0 (p) T q;ij (p, q) S 0 (q) where is the free fermion propagator, and is the effective vertex corresponding to the multiple scattering of the quark (or antiquark) off the shock wave background field. The dependence on the latter is given by the Wilson line defined in Eq. (3); here i and j label the colors of the incoming and outgoing quarks. Because we are including the possibility of "no scattering" within the definition of the effective vertex, the dressed propagator in Eq. 24 also contains a free part given by (2π) 4 δ (4) (p − q)S 0 (p) and an interacting part which contains all possible scattering with the nuclear shock wave. The diagram labeled as LO:(1) can be written as M LO: (1) µα;ij (q, k, p, k γ ) = −(eq f ) 2 d 4 l (2π) 4 u(k)γ α S 0 (k + k γ )T q;ik (k + k γ , q + l)S 0 (q + l)γ µ S 0 (l)T q;kj (−p, l)v(p) , (27) The integration over l − is trivial because of the δ(l − + p − ) factor from one of the effective vertices. The integration over l + is performed using the theorem of residues. After subtracting the "no scattering 1" contribution in which neither the quark or antiquark cross the nuclear shock wave, we can write the result compactly 5 as with 5 In the following, we will use the shorthand notations: The other diagram in Fig. 8 can be expressed similarly with R where v LO: The remaining R-functions are related to those in Eqs. 29 and 30 by the following replacements: If we keep the internal momentum labels identical to that in Fig. 8, this also results in an overall change in sign. Finally, one needs to redefine l ⊥ → −l ⊥ + k γ⊥ in order to make the transverse phases in all four contributions identical. The sum of the four contributions to the LO amplitude can therefore be compactly written as where q f is the charge of a quark or antiquark of a certain flavor f and is the sum of the contributions from the four processes, whose individual contributions are given by R LO:(i) µα . Plugging this expression for the amplitude (and its complex conjugate) back into Eq. (22), we obtain the LO triple differential inclusive cross-section for the production of a prompt photon in association with a dijet as [44] where α em = e 2 /4π is the electromagnetic fine structure constant, y = q · P N /l · P N is the familiar inelasticity variable of DIS and L µν is the lepton tensor in Eq. 21. We also introduced the differential phase space measures, d 6 K ⊥ = d 2 k ⊥ d 2 p ⊥ d 2 k γ⊥ and d 3 η K = dη k dη p dη kγ . In deriving the triple differential cross-section, we also isolated the prefactors (eq f ) 4 N c of the hadron tensor in Eq. 22 and used a properly normalized wave packet description for the incoming virtual photon [39,44]. The leading order hadron tensor is given bỹ where we introduced a compact notation for the integrals over the phases appearing in the amplitude expression in Eq. 32 The second such term appearing in Eq. 35 results from the complex conjugate of Eq. 36 and corresponds to replacing all transverse coordinates and internal momenta therein by their primed counterparts. The function represents the spinor trace in the cross-section. Finally, the nonperturbative input from the dynamics of saturated gluons in the nuclear target is contained in Ξ(x ⊥ , y ⊥ ; y ⊥ , x ⊥ ) which can be decomposed as where represent respectively dipole and quadrupole Wilson line correlators. A pictorial representation of these correlators is given in Fig. 9. These gauge invariant quantities appear in a variety of processes in both p + A and e + A collisions. Explicit expressions for these correlators are available [8,89,90] in the McLerran-Venugopalan model [3][4][5], where the distribution of sources Here where µ 2 A = A/2πR 2 ∼ A 1/3 is the average color charge squared of the valence quarks per color and per unit transverse area of a nucleus with mass number A. As we will discuss later, the JIMWLK evolution equation can be reexpressed as evolution equations for these gauge invariant quantities.

III. OUTLINE OF THE NLO COMPUTATION
Before we dive into the rather involved computations (which, as articulated briefly in the introduction, have much of the complexity of two-loop computations in standard pQCD) it is useful to outline the structure of various contributions to the computation at NLO. These can be classified into real and virtual contributions; the latter can be further subdivided into self-energy and vertex corrections. An important simplification in the Regge limit is that the shock wave interaction is instantaneous, which eliminates more than one insertion from the effective vertex on any given line in a Feynman diagram. In addition to outlining the structures comprising the different contributions, we will also provide in this section a flow chart which points to the different contributions, and links that take the reader to specific terms in the computation, without having to wade through the entire detailed computation in the next section.
A. Structure of contributing processes There are both real gluon radiation and virtual gluon exchange processes that contribute at O(α S ) to inclusive photon+dijet production. For the computation of the NLO differential cross-section, we need to take the modulus squared of the amplitudes for gluon emission for fixed static color sources, perform the CGC averaging over the distribution of these color charge configurations, and finally integrate over the phase space of the emitted gluon. In the case of the amplitudes of graphs containing virtual gluons, we need to include the interference of these with the leading order amplitude given by Eq. 32, before performing the CGC average over static color charge distributions.
The NLO hadron tensor is then given by where is shorthand notation for integration over the phase space of the emitted gluon and c.c denotes the complex conjugate. For virtual exchange graphs, we can broadly classify the two topologies of diagrams as self-energy and vertex contributions, which we have denoted above with the superscripts 'SE' and 'Vert' respectively. We will now describe the further systematic classification of the contributions to the amplitudes in each category in terms of their color structure.

Real emissions
There are 20 Feynman graphs that describe the radiation of a gluon in addition to the photon radiated in the final state. Further, there are distinct topologies of these graphs depending on whether 1. the gluon is emitted prior to scattering of the quark and antiquark, or 2. emitted by the quark or antiquark after they scatter off the nucleus.
In the former case, the gluon has the possibility of scattering off the background classical field whereas in the latter case it does not. For each of these diagrams, we need to subtract the "no-scattering" contribution to the amplitude, which is obtained by setting U andŨ 's to unity.
As in the case of the LO amplitude, we can write the amplitude for real emissions as where dΠ R represents the integrals over the transverse Fourier phases associated with the effective vertices. By an appropriate redefinition of momenta, these can be made identical for all the contributions. Their exact form is not important for the present discussion but will be delineated in the upcoming sections which contain the detailed computation of the amplitudes for the various processes. The essential features of T R , and T R are as follows: • There are a set of 10 diagrams that contribute to the factor T R . These are the processes where the emitted gluon may simultaneously scatter off the background classical field in addition to scattering of the quarkantiquark dipole. A representative diagram is shown in Fig. 10, with the other diagrams obtained simply both by permutations of the emission vertex for the final state photon and by interchanging the quark-antiquark lines. • There are 5 contributions that constitute T (2) R in Eq. 44. These correspond to gluon emission from the quark after it scatters off the nucleus. A representative graph is shown in Fig. 11; the others are obtained by permutations of the vertex for the final state photon.
• Finally, there are 5 contributions constituting T (3) R which are obtained by interchanging the quark and antiquark lines in the Feynman graphs of Fig. 11. These are identified separately because they have a different color structure from the diagrams comprising T Thus at the level of the squared amplitude, 400 diagrams contribute to the NLO photon+dijet cross-section.

Virtual contributions
Broadly speaking, virtual contributions can be classified into vertex and self-energy graphs. In addition, there are diagrams in which the emitted gluon scatters off the shock wave before being reabsorbed by the quark/antiquark. To add to the complexity of such computations, the photon can be emitted either before or after these scatterings from the quark or antiquark. Thus the total number of diagrams to compute is significantly more than fully inclusive DIS at NLO. These can however be classed into distinct categories based on their Wilson lines structures. 1. Self-energy contributions: We can write the amplitude of the self-energy contributions as where dΠ S denotes the phase space factor corresponding to the self-energy contribution.
• There are 6 contributions proportional to T (1) S in the expression above; one such diagram is shown in Fig. 12. The topology of these diagrams corresponds to that of a gluon emitted by the quark prior to scattering and then reabsorbed by the quark after the qqg state scatters off the shock wave.
• The 6 contributions that constitute T  • There are 24 other contributions in which the emission and absorption of the gluon occurs either prior to or subsequent to scattering with the nucleus; these are proportional to T  2. Vertex contributions: Similarly to the self-energy contributions, the general expression for the amplitude of vertex contributions can be written as where dΠ V represents the phase space factor for vertex-like corrections and C F = (N 2 c −1)/2N c is the quadratic Casimir for the fundamental representation of SU (N c ).
• There are 6 contributions to T (1) V . A typical diagram is shown in Fig. 14; the rest are obtained by permutations of the photon emission vertex. These correspond to the virtual gluon emitted by the antiquark following which it crosses the shock wave before being absorbed by the quark.
• The T • There are 6 contributions proportional to T V ; one such graph is shown in Fig. 15. These are part of the radiative corrections to the virtual photon wavefunction fluctuating into a quark-antiquark dipole with the addition of a final state photon. Consequently, the Wilson line factor is identical to that in the LO amplitude times the color factor C F .  • Finally, we have 6 contributions proportional to T  For the convenience of the reader, the computational tree depicted in Figure 17 shows the components and subcomponents building up the NLO hadron tensor in Eq. 42. Clicking on each of these will take him or her to the particular expression desired. This computational tree is therefore also a template for numerical evaluation of the NLO photon+dijet cross-section that will be the subject of future work.
As discussed above, perturbative contributions from kinematically allowed diagrams with similar color structure are contained in the various T R,S,V functions. Each of these is the sum of the contributions of the different Feynman diagrams denoted by R i µα and are presented in columns under the T R,S,V functions in Fig. 17. Within a certain column, there may be diagrams that are connected to one another by quark-antiquark interchange. We have put these together within blue rectangular boxes in Fig. 17 . These R-functions can be obtained in sequence by imposing the q ↔q replacements given by Eq. 66 (later in the text) in the functions appearing in the columns above the blue boxes. Moreover, there are entire categories of processes related by interchange of the quark and antiquark lines. These are also shown in Figure 17.

B. Assembling the different contributions in the amplitude squared
For the computation of the differential cross-section, we need to take the modulus squared of the amplitude for the real emission processes and the interference of the virtual graphs with LO processes. The general expressions for the NLO amplitudes are given by Eqs. 44, 45 and 46 respectively while Eq. 32 denotes the same for the LO amplitude.
FIG. 17. Computational tree for the NLO computation of the hadron tensor. The different branches correspond to the different components constituting X NLO µν . The first of these nodes represent the amplitude contributions that comprise the NLO hadron tensor. The sub-branches contain the combined result for the hard parts of the amplitude for the different contributing processes. These individual contributions from diagrams which are categorized based on the color structure are provided in the long columns. Quark-antiquark interchanged counterparts of quantities appearing in the same class of diagrams and hence the same column are grouped within blue rectangular boxes. As an example, R (R6),...,(R10) are obtained respectively by exchanging q ↔q in R (R1),...,(R5) . The labeling follows the same ordering for the other quantities grouped in blue boxes. Terms grouped in the red box are all zero and do not contribute to the amplitude.
The squared amplitude, a functional of the stochastic source charge density ρ A , then needs to be averaged over all possible charge configurations weighted by the distribution W [ρ A ]. Following extensive use of the Fierz identity and the relationŨ connecting adjoint Wilson lines to fundamental ones, we get non-trivial combinations of dipole and quadrupole Wilson line correlators (see Eqs. 39). These are summarized clearly in Table II below. The terms proportional to products of Wilson line factor Real emission Virtual: Vertex Virtual: Self-energy T's represent Dirac traces obtained from expressions for the squared amplitudes. The corresponding color trace over products of Wilson lines corresponding to each row is given in the leftmost column. To obtain the color factors for the conjugates of the terms in rows 4 (5), we need to replace (x ⊥ , y ⊥ ) → (y ⊥ , x ⊥ ) in the transverse coordinates of the corresponding color factors of rows 5 (4). As evident from Table II, the fundamental building blocks which span the entire high energy computation have the structures D, Q, DD and DQ, albeit with different dependence on the transverse coordinates. In the sections that follow, we will carry out detailed computations of the various entities in Table II. The organization of the NLO computation in the manner described here will provide a transparent guide to the identification of soft, collinear and ultraviolet divergences in the computations.

IV. NLO CONTRIBUTIONS TO THE AMPLITUDE FROM REAL EMISSIONS: DETAILED CALCULATIONS
In this section, we will compute in detail the amplitudes for the various real emission graphs presented in Sec. III. As discussed there, there are two distinct topologies based on gluon emission before and after scattering of the dipole off the background classical field. We shall now systematically illustrate how to calculate the various diagrams contributing to each of the three terms in the general amplitude expression in Eq. 44. Readers uninterested in these details can proceed directly to Section VI.

Contributions to T
(1) R : The processes that contribute to T (1) R in the amplitude (see Eq. 44) for real gluon emission are shown below in Fig. 18. Only half of them (labeled (R1) − (R5)) are presented here. The quark↔antiquark interchanged counterparts of these 5 diagrams are respectively labeled (R6) − (R10). Before delving into the details of the computation, we will write down the general form of the contribution from these 10 processes to the total amplitude. After subtracting the "no scattering" contribution from each of them, this is given by 18. Real emission diagrams contributing at NLO to the "impact factor" with the gluon emitted prior to the scattering of the quark and antiquark off the nucleus. The other five diagrams are labeled (R6) − (R10) and can be obtained respectively by quark-antiquark interchange of (R1) − (R5). These 10 contributions constitute in total the coefficient T Here z r tot = z q + zq + z γ + z g is the total momentum fraction for real emission and dΠ LO ⊥ is a shorthand for the integrals We will now discuss the computation of the R (Rβ) 's that constitute T (1) R given by Note that in this discussion, and all subsequent discussions, we will only explicitly show the dependence (if any) of these functions on internal momenta (that are integrated over) albeit they are of course functions of the external momenta as well.
The contribution to the amplitude for the processes labeled (R1) in Fig. 19 is given by 19. The NLO process labeled (R1) in Fig. 18 with all momentum assignments and directions shown. The effective vertices are clearly shown in Fig. 1.
where the free fermion and gluon propagators in A − = 0 gauge are respectively given by The effective vertices for the quark and gluon are contained in the expressions for their dressed propagators in the background of the strong classical color field of the nucleus. Recall that the expression for the quark propagator is (given previously in Eq. 24) In the "wrong" light cone gauge A − = 0, we conveniently obtain an analogous form for the dressed gluon propagator which can be written as [5,15,[68][69][70][71]] where µ, ν and a, b are the Lorentz and adjoint color indices for the outgoing and incoming gluon which respectively carry momenta p and q. Again, to recapitulate, the expressions for the effective vertices (introduced in Fig. 1) are, Since n. * (k g ) = 0 in A − = 0 gauge, the polarization vector for the outgoing gluon has the form (k g ) = Using this, we can derive the following useful relation where we have used the eikonal approximation l − 2 = k − g contained in the expression for the effective gluon vertex. Now integrating over l − 1 and l − 2 using the δ-functions appearing in the effective vertex factors, we can rewrite Eq. 52 as where r zx = z ⊥ − x ⊥ and the numerator and denominator are respectively given by and As expected, we have an overall longitudinal momentum conserving δ-function where P − tot = k − + p − + k − γ + k − g , which is a reflection of the eikonal approximation ingrained in our analyses. When we examine the above equations, it is clear that the numerator structure allows for the use of Cauchy's residue theorem to evaluate the + -integrals by complex contour integration. There are two l + 2 poles on either side of the real axis. We deform the contour clockwise so as to enclose the pole below and subsequently perform the l + 1 integration by an anticlockwise contour deformation. We next perform the l 2⊥ integration using the expressions for the relevant integrals tabulated in Appendix B (see Eqs. B9 and B10 for the expressions in d dimensions). Finally, subtracting the "no scattering" contribution by setting the Wilson lines to unity, we write the resulting amplitude as where is independent of k g⊥ but depends on the other external momenta which are not shown explicitly in its argument.
In the above equation, the functions I are proportional to modified Bessel functions of the second kind (or Macdonald functions). In d = 2 dimensions, and for the process (R1), these can be written as with the arguments of the functions given by At the level of the differential cross-section, we will integrate over the phase space of the real gluon which includes an integration over k − g from 0 to +∞. If we examine closely Eqs. 63 and 64 above, we observe a logarithmic singularity in the limit k − g (z g ) → 0. The other limit (k − g → +∞) converges because of the oscillatory nature of the exponentials in the I-functions. We will show later that this slow gluon limit (z g → 0) is what generates the large logs in x -the net contribution of terms with these large logs multiplies the JIMWLK kernel. This aspect of the computation will be discussed at length in Sec. VII. The remaining four diagrams in Fig. 18 can also be computed in a similar fashion and the combined contribution is given by where the R-functions are given in Appendix C. In order to find the corresponding contributions of Fig. 18 (with the quark and antiquark lines interchanged) which we call (R6) − (R10), we need to impose the following replacements in the R-functions in Eqs. 62 and C2-C9 As discussed at the end of Sec. II, the last redefinition is only to ensure that the transverse phases defined by Eq. 50 remain the same so that the net contribution to the amplitude from the 10 diagrams can be compactly written as in Eq. 49.

Contributions to T
(2) FIG. 20. Real emission graphs with the gluon emitted by the quark subsequent to its scattering off the nucleus. The graphs obtained by interchanging the quark and antiquark lines in the above diagrams are respectively labeled (R16) − (R20) and they constitute T We will now compute the contributions from diagrams shown in Fig. 20 in which the gluon is emitted by the quark after it scatters off the background classical field. The combined amplitude from these 5 processes can be written as where has an explicit k g⊥ dependence.
In the following, we will show how to compute one such contribution. The rest can be computed using the same techniques. The contribution to the amplitude from the diagram labeled (R11) , with the detailed momentum assignments and directions shown in Fig. 21, is given by where l − 1 has been integrated out using the δ-functions contained in the expressions for the effective vertices in Eqs. 56. We perform the integration over l + 1 by a clockwise deformation of the contour. Finally subtracting the "no scattering" contribution, we get the amplitude as Fig. 20 with the momentum assignments and directions shown. The gluon does not suffer scattering off the nucleus in the above scenario. where At the level of the inclusive cross-section when we integrate over k g⊥ , it is evident from Eq. 71 that we will once again encounter tensor integrals of the kind given by Eq. B8. The other contributions arising from the emission of the gluon after the scattering off the shock wave can be similarly computed and their combined contribution is represented by Eq. 67. Expressions for the remaining R-functions are provided in Appendix C.

Contributions to T
R : These are the processes obtained by interchanging the quark and antiquark lines in (R11)−(R15). They are respectively labeled (R16)−(R20) and their contributions to the amplitude are obtained by imposing the replacements given by Eq. 66 in Eq. 68. This ensures that the transverse phases remain the same throughout. The total contribution from this final sub-category of diagrams can then be written as with Finally, we can write the combined contribution from all the allowed 20 real emission diagrams as To compute the contributions of these real graphs to the differential cross-section, we need to take the modulus squared of the amplitude in Eq. 75 and then integrate over the phase space of the real gluon. From the discussion here, and the expressions given in Appendix C, it is clear that for the 10 diagrams (see Fig. 18) in which the gluon gets scattered off the nucleus, the amplitudes can be written in terms of MacDonald functions whose arguments in general depend on the gluon momenta. As such, it is difficult to isolate analytically the rapidity divergent pieces and the finite contributions from the squared amplitudes for these graphs. In Sec. VII, we will obtain the rapidity divergent pieces from these amplitudes by explicitly taking the k − g → 0 limit and show that these pieces contribute towards small x evolution. To compute the finite pieces from this class of diagrams, one however needs to perform the integration over the gluon phase space numerically by imposing a cutoff for the gluon momentum fraction z g . Because of the interaction with the nuclear shock wave, there are no collinear divergences associated with these diagrams.
For the processes shown in Fig. 20 (and their q ↔q counterparts) in which the gluon does not scatter off the nucleus, there are divergences from the region of phase space where the gluon is soft (k g → 0) and/or collinear (k g⊥ ∝ k ⊥ , p ⊥ ) respectively to the antiquark and quark. In Sec. VI, we will promote partons to jets and explicitly extract these divergent structures by using a jet cone algorithm. This will allow us to show the cancellation between residual collinear divergences from the virtual graphs with those in the real gluon amplitude squared and therefore obtain an IR safe cross-section.

V. NLO CONTRIBUTIONS TO THE AMPLITUDE FROM VIRTUAL GRAPHS: DETAILED CALCULATIONS
In this section, we will illustrate the details of the computation of the amplitudes corresponding to the virtual diagrams shown in Sec. III. We will start with the self-energy diagrams and follow this with the computation of the vertex correction graphs. An additional feature of these processes relative to the usual Feynman diagram computations is that the emitted gluon can scatter off the background shock wave classical field before being absorbed by the quark or antiquark.
A. Self-energy graphs with dressed gluon propagator As discussed previously, there are three distinct topologies of the Feynman graphs describing self-energy contributions. These are discussed individually below.

Contributions to T
(1) S : The diagrams contributing to T (1) S in the general expression for the amplitude given by Eq. 45 are presented in Fig. 22. These are the processes which allow for a virtual gluon emitted from the quark line to scatter off the shock wave before being reabsorbed. We will first present the combined result for the amplitude from all such processes and then demonstrate the details of the computation using a representative diagram.
The combined amplitude from these 6 processes has the structure where z v tot = z q + zq + z γ is the total momentum fraction for virtual emission and we have bundled together the transverse Fourier phases using the short-hand convention which contains an additional integration over z ⊥ relative to the similar expression given in Eq. 50. Finally,

T (1)
S;µα can be written as the sum of a piece that includes UV and rapidity divergences and a finite part, Once again, we are including the dependence on momenta that are integrated over. Also note that as previously, r zx = z ⊥ − x ⊥ . In the above equation, T LO (given previously in Eq. 33) contains contributions to the amplitude from the four allowed LO processes. The pieces that are independent of l 2⊥ will produce a δ (2) (r zx ) from the l 2⊥ integration. Once the z ⊥ integration is done, the color structure corresponding to these pieces then reduce to that for the LO amplitudes times the quadratic Casimir C F . We can therefore write the amplitude in Eq. 76 as where M LO is the leading order amplitude given by Eq. 32.
The divergence free pieces of the various contributions are combined to constitute M SE (1) finite . The logarithmic divergence ln(1/z 0 ) = ln(q − /Λ − 0 ) in the above expression, where Λ − 0 is given by Eq. 4, arises from the integration over the momentum fraction z l of the gluon in the loop. Its upper limit, up to logarithmic accuracy at this order, is controlled by the q − momentum component of the photon and its lower limit by the longitudinal width 1/Λ − 0 (or equivalently 1/P + N ) of the target nucleus. The latter is a cutoff that we are imposing to regulate the gauge pole l − = 0 in the LC gauge gluon propagator.
Further, from the expression for the amplitude in Eq. 79, we can see that there are two kinds of singular logarithms multiplying ln(1/z 0 ). The one that appears in the first line of Eq. 79 arises from the collinear limit z ⊥ → x ⊥ and are not part of the small x logarithms contributing to JIMWLK evolution. Some of these divergences will cancel out already at the amplitude level between different classes of diagrams and the rest will cancel between real and virtual graphs. The limit z ⊥ x ⊥ , y ⊥ of the evolution kernels is captured by the logarithms appearing in the second line of Eq. 79. At this order of accuracy, the ln(1/z 0 ) log can also be expressed as the ln(x 0 /x Bj ) log giving rise to small x JIMWLK evolution. We will discuss this point in greater detail in Section VII. In the limit of large Q 2 , these will give the double log limit of the DGLAP evolution equation [91][92][93][94] or equivalently the large Q 2 limit of BFKL equation. The 1/ singularities for → 0 arise from regulating the UV divergences in the integrations over the transverse loop momentum of the gluon using dimensional regularization in d = 2 − dimensions. In our expressions, We will now use one representative process to explicitly demonstrate the steps leading to the above result, in particular, the computation of the divergent pieces. The calculation of the finite pieces, while absolutely essential for precision computations, are not particularly illuminating; they are discussed in Appendix I. The amplitude for the process labeled (S1) with momentum assignments shown in Fig. 23 is given by 23. The process labeled (S1) in Fig. 22 with momenta and their directions shown. We must subtract the "no scattering" case (corresponding to the same amplitude without crossed or filled blobs) from the above contribution to obtain the physical amplitude.
where the free fermion and gluon propagators are given respectively in Eq. 53 and the corresponding effective vertices are in Eqs. 56.
Although the structure of the free gluon propagator in A − = 0 gauge is in general more complicated compared to that in covariant gauges, we will use the light cone condition to derive an identity that allows one to simplify the amplitude above. To show this, we start with the expression where (. . .) denotes the terms between the two gamma matrices. The terms in parentheses on the right denote the Lorentz structure of the free gluon propagators while the metric g σρ in between these is from the effective gluon vertex. Using n µ = (1, 0, 0 ⊥ ) and some algebra with the indices, it is possible to reduce the above expression to where k = 1, 2 is summed over. Using the above identity, analogous to the one in Eq. 57, and integrating out l − 1 , l − 2 using the δ-functions in the effective vertices, Eq. 80 can be written as with and .
Note that in these expressions we have canceled a factor of 2l − 3 from the gluon effective vertex in the numerator with a corresponding factor from one of the propagators in the denominator. From the above expression for the denominator, it is evident that there are two l + 3 poles which are located at implies that the poles are on the same side (above) of the real axis. Since the numerator is independent of l + 3 , the contour can always be deformed such that none of the poles are enclosed, thereby giving a null result. Hence we must have l − 3 > 0 as well as k − + k − γ − l − 3 > 0 for a non-zero contribution. As discussed earlier, logarithmically divergent integrals in l − 3 are regulated by introducing a lower cutoff Λ − 0 given by Eq. 4. We observe that there is an equivalence of this picture to analyses in light cone perturbation theory, where the positivity of l − 3 here correspond to forward propagation (in light cone time) of the exchanged gluon. From inspection, an identical argument holds for the l + 2 pole. We will enclose l + 3 | B and the following poles for the contour integration over l + 1 and l + 2 respectively, Finally subtracting the "no scattering" contribution we arrive at the following expression for the amplitude, where dΠ v ⊥ is defined by Eq. 77, z l = l − 3 /q − is the gluon momentum fraction in the loop and v (S1) As previously for the real contributions, we have here too the familiar integrals over transverse momenta l 2⊥ and l 3⊥ , that we encountered in the computation of the real emission amplitude, which contain the exponential terms exp(±iv (S1) i⊥ .r zx ) (i = 1, 2) and are proportional to Macdonald functions. We also learnt from our previous discussion that these can give singular logarithms from the z l integration in the limit z l → 0 when these phases become unity. We will see later that this is indeed the case. Another interesting possibility is the limit r zx → 0 which corresponds to the UV limit for the momenta, l i⊥ (i = 2, 3) conjugate to this transverse coordinate vector. The manifestation of this UV divergence is explicit if we perform the following momentum redefinition, l 2⊥ − l 3⊥ → l 2⊥ in Eq. 88 with l 3⊥ remaining unchanged.
If we now do a naive power counting in l 3⊥ , it is clear that the numerator (N (S1) ) and denominator (D (S1) ) are both proportional to l 4 3⊥ for large l 3⊥ . By a careful use of Dirac algebra, it can be shown that the terms proportional to l 4 3i and l 3 3i vanish and hence the transverse momentum integral is at most logarithmically divergent. We will use dimensional regularization 6 in d = 2− dimensions for terms in the integrand proportional to l 2 3⊥ and use the d = 2 results for the convergent pieces. However we also have to account for rapidity logarithms from the z l integration which as one can see arises from terms proportional to 1/z l . We will therefore separate our amplitude into two contributions. The first contribution is constituted of terms in the integrand proportional to l 2 3⊥ and also includes terms which will be finite from the l 3⊥ integration but contains 1/z l pieces . The second contribution is comprised of terms proportional to l p 3⊥ (p < 2) and does not contain any piece proportional to 1/z l . This will then be a genuinely finite part of the amplitude.
With this in mind, we can rewrite the amplitude in Eq. 88 as where R (S1) contains the UV and rapidity divergent pieces and some finite terms. We will isolate these remaining finite pieces and combine them with the genuinely finite contribution R (S1) (II) to obtain the net finite contribution from (S1). Carefully isolating these pieces from the numerator in Eq. 88 and after an extensive use of Dirac algebra in d-dimensions and the identity γ i γ j l 3i l 3j = −l 2 3⊥ we can simplify the first contribution to read R (S1) where the constituent integrals appearing above are given by and We have also introduced here the short-hand notations where the relevant terms v (S1) 1,2⊥ and ∆ (S1) 1,2 were previously defined in Eqs. 89. Using the standard identities for → 0 in Eq. 91 and definingμ 2 = 4πµ 2 e −γ E , we can finally write where we have deliberately separated the logarithm containing l 2⊥ by introducing the resolution scale Q 2 . We will show later that divergent terms of the kind appearing in the first line of the above equation cancel and do not appear in the final cross-section. This in turn also shows the independence of our final observable on the scaleμ.
The finite pieces from R (S1) (I) that do not contain any logarithms in z 0 nor UV logs are contained in the remainder term defined as 7 (S1) In Appendix I, we show the computation of the above term in detail. The R-function appearing in Eq. 90 can finally be written as In case of the diagrams (S2) and (S3), we get two independent contributions to the amplitude because we have different choices of contours for the integration over l + 2,3 depending on whether 0 < l − The net amplitude is therefore obtained by summing these individual contributions. We can show that for either range of l − 3 the divergent pieces in R (Sβ) (I);µα (β = 2, 3), are proportional to g α+ which will yield zero when contracted with the polarization vector for the outgoing photon because − * (k γ ) = 0 in our choice of gauge. Therefore, for these diagrams, only the finite pieces survive; we will present the expressions for these in Appendix I.
Finally, for diagrams (S4) − (S6), we can show that the divergent parts of the amplitudes have exactly the same structure as Eq. 96 albeit with the (previously specified) different R LO 's depending on the topology of the diagram.
We would like to remind the reader that in our chosen convention the q ↔q interchanged counterparts of the LO processes labeled LO:(1) and LO:(2) as shown in Fig. 8 are respectively labeled LO:(3) and LO:(4). We can now finally express the sum of the contributions to the amplitude from the six processes (in Fig. 22) in the form given by Eq. 79.
The finite pieces of the amplitude are contained in where These are presented in detail in Appendix I.
The fact that the divergent part of the amplitude in Eq. 79 is proportional to the LO amplitude shows that scattering off of the background classical field does not affect the UV structure of these processes. This is to be expected because the gluon in the loop only experiences transverse momentum kicks while propagating through the nuclear shock wave; in the limit of large loop momentum, this scattering should have no effect on the short distance structure of the theory.

Contributions to T
(2) S : These are the processes labeled (S7) − (S12) which are obtained by interchanging the quark and antiquark lines in the diagrams shown in Fig. 22. The combined contribution to the amplitude from these processes can be written as where T (2) S is obtained by imposing the replacements in Eq. 66 to the corresponding expression in Eq. 78 of T (1) S . The resulting expression is given by The amplitude in Eq. 104 can therefore be rewritten as where the finite piece is with As noted earlier, the two functions appearing above are obtained from their quark↔antiquark interchanged counterparts by imposing the replacements in Eq. 66.

B. Self-energy graphs with free gluon propagator
Contributions to T S : There are a total of 24 diagrams which contribute to the quark self-energy corrections at NLO. Half of them are shown in Fig. 24 and the other half are obtained simply by interchanging the quark and antiquark lines. We have grouped them in three rows depending on their divergence structure. As we will show, the first four processes labeled (S13) − (S16) inherit a UV divergence structure which is identical to the one described in the previous section.
Because we are working in the limit of massless quarks, we do not need to renormalize the quark mass. So diagrams in which the gluon loop is on an external on-shell quark or antiquark line will contribute zero to the amplitude. This is true for the four diagrams labeled (S17) − (S20) appearing in the second row of Fig. 24. A detailed computation of the quark self-energy in Appendix D in the "wrong" LC gauge A − = 0 demonstrates this result.
We will also use the expressions obtained in this Appendix for the self-energy loop to compute the diagrams which have a similar topology albeit different locations of the emission of the photon. Finally we will compute here the amplitude for the two processes (S22) and (S24) which represent α S corrections to the quark-photon-quark vertex. Although these diagrams will yield a different UV divergence structure from the rest, we will see that there are cancellations amongst the divergences of the four processes in the third row (which explains why they are grouped together). A general expression for the gluon loop correction to the quark-photon-quark vertex is derived in Appendix E.
We will begin our discussion by considering the amplitudes for (S13) − (S16). Subtracting the "no scattering" contribution from each of these diagrams, we can write the sum of the amplitudes from these processes as In the above equation, T LO is given by Eq. 33 and the coefficients A (Sβ) and B (Sβ) will be given later. The finite pieces of each process are contained in R (Sβ) finite and expressions for these are provided in Appendix I 3.
We shall now discuss in detail the computation for (S13). The results for the other three can be obtained using the same methods. The amplitude for this diagram (shown in Fig. 25) is 25. The process labeled (S13) in Fig. 24 with the momentum assignments and directions shown. The gluon does not scatter off of the background classical field in this picture. Only the quark and antiquark do which are represented by the crossed blobs.
M (S13) where the free quark propagator and effective vertices are given respectively in Eqs. 53 and 56. The self-energy contribution Σ computed in Eq. D10 of Appendix E is with the four-momentum k f = (q + + l + 1 , k − + k − γ , l 1⊥ ). The constituent integrals are given by We have used here the conservation of momentum across the eikonal quark vertex to write q − + l − 1 = k − + k − γ . The integration over l − 1 follows trivially from the Dirac delta functions contributed by the effective vertices and we will obtain an overall momentum conserving delta function. Since k 2 f depends on l + 1 the value of ∆ s will be determined by the pole enclosed in the contour integration over l + 1 . We will enclose the pole at which is above the real l + 1 axis. The ∆ s factor, appearing in the constituent integrals in Eqs. D11, for (S13) therefore becomes ∆ (S13) where z l = l − 2 /q − is the gluon momentum fraction in the loop and ∆ LO:(1) = Q 2 zq(1 − zq) − iε. Using the identities given by Eqs. D14 and D11 in Eq. 111, we can finally write the expression for the expression in Eq. 111 as where we employ the compact notation for the transverse Fourier phases appearing in the LO amplitude introduced in Eq. 50. For clarity, as previously, we split R (S13) as where the first part is constituted of logarithms in z 0 (equivalent to logarithms in rapidity) and UV singularities and the second part is free of such terms. For (S13) these are obtained respectively to be and R (S13) Following the structure of the amplitude for these four processes given by Eqs. 109 and 110, the coefficients A and B for (S13) are respectively We can now employ these techniques to compute the contributions (S14) − (S16): The A and B coefficients above are where and Finally we have The terms v LO:(1,2) ⊥ and ∆ LO: (1,2) are given in Eq. 31 and ∆ LO: The finite pieces for these processes are given in Appendix I.
Moving now to the third line of Fig. 24, we can use the self-energy function in Eq. D10 to compute the contributions to the amplitude from diagrams (S21) and (S23). As for the other diagrams, the expression for the amplitude can be written as where We will provide expressions only for the divergent parts of these amplitudes in this section. The finite pieces are given in the Appendix I 3. For (S21) and (S23) the divergent pieces are respectively given by and We conclude this section by detailing the computation of the amplitude for the remaining two diagrams labeled (S22) and (S24) that appear in the third row of Fig. 24. For this we will use the general results in Appendix E for the gluon loop contribution. For the diagram (S22), we write the amplitude as In Appendix E, we wroteΣ α as the sum of two contributionsΣ =Σ A α +Σ B α where A and B denote respectively the cases where the gluon loop momentum l − is in the range (Λ − 0 , k − ) and (k − , k − + k − γ ). The divergent pieces of each of these contributions are given respectively in Eqs. E10 and E24. Using the Dirac equation u(k)/ k = 0 for massless quarks, and a bit of algebra involving gamma matrices, we can rewrite Eq. 130 as where The divergent piece is We can clearly see from Eqs. 128 and 133 that the divergent pieces for (S21) and (S22) exactly cancel leaving only the finite pieces from these graphs towards their contribution to the amplitude. In a similar fashion, we write the amplitude for (S24) as where The divergent piece for this diagram is provided in Eq. E25 of Appendix E. Comparing this expression with the divergent part of (S23) (see Eq. 129), we can check that there are indeed UV and rapidity divergence cancellations occurring in the net contribution from these processes. The only UV divergent term that survives is the one whose Dirac gamma matrix structure is γ − γ α γ µ . As we will show in the next section, a similarly divergent piece is obtained from the vertex correction graph labeled (V 15) in Fig. 30 that cancels this UV divergence. The result for the finite piece of the loop diagram in Fig. 34 is given in Appendix E. This result has been used to compute the finite pieces of the amplitude for (S22) and (S24) in Appendix I 3.
We can now write down the final expression for the amplitude from the 12 processes shown in Fig. 24 (and their quark↔antiquark interchanged counterparts) that contribute to the "free gluon" self-energy amplitude as where The (q ↔q) interchange in the above equation represents the contribution to the amplitude from the processes labeled (S25)−(S36). These are obtained by imposing the replacements in Eq. 66 in the result we computed for (S13)−(S24). The divergent part of the amplitude coming from the 12 diagrams in Fig. 24 is given by (2) . (138) The terms appearing in the first line of the above equation are equal in magnitude but opposite in sign with respect to the amplitudes from the six diagrams (S1) − (S6) given by Eq. 78. They will therefore cancel out when we add the contributions from all possible self-energy corrections.
There are two kinds of rapidity divergent pieces in the second line of the above equation. (We have used µα (l 1⊥ ) = T LO µα (l 1⊥ ) in writing the second term.) The terms proportional to ln(1/z 0 ), with the coefficients A (Sβ) given in Eqs. 119, 121, 122 and 125 contribute to the LO JIMWLK Hamiltonian. Indeed, these are precisely the terms that will give the double log limit of the DGLAP/BFKL equations in the limit of large Q 2 . On the other hand, the terms proportional to ln 2 (1/z 0 ) cancel between real and virtual graphs; we will demonstrate this explicitly in section VI.
Finally, the term in the last line of Eq. 138 is the UV divergent piece that remains after cancellations between the divergent pieces of graphs (S23) and (S24). This divergence will be canceled by a similar contribution from the vertex correction graph (V 15) which also has a photon nested in the gluon loop.

C. Vertex graphs with dressed gluon propagator
In this section, we will compute the vertex corrections to the LO amplitude in which the gluon crosses the nuclear shock wave. As discussed in Sec. III, there are two distinct topologies with six contributions in each class. These two sets of contributions are related to each other by q ↔q interchange. Following a logic identical to the discussion of self-energy graphs, we will detail the computation of the amplitude for one such representative process. The remaining processes can then be computed following similar techniques.

Contributions to T
(1) in Eq. 46. The topology corresponds to that of a gluon emitted by the antiquark which propagates forward and is absorbed by the quark after the qqg state scatters off the nucleus. The different diagrams correspond then to the possible locations from where the final state photon might be emitted. By interchanging the quark and antiquark lines, we get the diagrams constituting T The six diagrams contributing to T (1) V in the general structure of the amplitude given in Eq. 46 are shown in Fig. 26. Their combined contribution to the amplitude can be written as where z v tot = z q + zq + z γ and is the net perturbative contribution from the six processes which can be expressed as the combination of a divergent and a finite part. We will see that there are no UV divergent pieces for these diagrams and the only singularities are those arising from the l − = 0 pole in the free gluon propagator. In the following, we will outline the steps leading to the above conclusion by considering the representative process (V 1).
The amplitude for (V 1) with the momentum assignments shown in Fig. 27 is given by 27. The process labeled (V 1) in Fig. 26 with momenta and their directions shown. We must subtract the "no scattering" case from the above contribution to obtain the physical amplitude. This will correspond to the diagram above but with no crossed and filled blobs.
where we remind the reader that the free fermion and gluon propagators are given respectively in Eq. 53 and the corresponding effective vertices in Fig. 1. The various Lorentz indices appearing above can be contracted to obtain the simplified expression given by Eq. 82 for the gamma matrices. The integrals over l − 1 and l − 2 can be trivially performed using the delta functions appearing in the effective vertices. The integrations over l + 1 , l + 2 and l + 3 can then be systematically performed using Cauchy's theorem of residues. These contour integrations are non-zero for l − 3 in the range 0 < l − 3 < (k − + k − γ ). Finally introducing the gluon loop momentum fraction, z l = l − 3 /q − and subtracting the "no scattering" contribution we can write the (V 1) amplitude as where the numerator and denominator are respectively given by and Note that in these expressions we have canceled a factor of 2z l from the gluon effective vertex in the numerator with a corresponding factor from one of the propagators in the denominator. This particular form of the amplitude given by Eqs. 142, 143 and 144 will be useful later in Sec. VII. For extracting the UV divergent pieces we will redefine l 2⊥ − l 3⊥ → l 2⊥ and l 3⊥ − z l /(1 − zq) (k ⊥ + k γ⊥ ) → l 3⊥ to rewrite Eq. 142 as where the modified numerator and denominator arẽ To write this denominator compactly, we have introduced the variables v (V 1) The denominator in Eq. 147 grows as l 6 3⊥ for very large l 3⊥ . As we argued in the computation of (S1) in Sec. V A, terms in the numerator in Eq. 143 proportional to l 5 3⊥ and above vanish. This implies that the integral over l 3⊥ is at most logarithmically divergent in the limit of large l 3⊥ . Collecting terms in the numerator proportional to l 4 3⊥ (and using γ i γ j l i 3 l j 3 = −l 2 3⊥ , γ − γ µ γ − = 2g µ+ = 2δ µ− ), we can write the piece of the numerator (in d = 2 − dimensions) contributing to the UV divergence as However from the general expression for the amplitude of our photon+dijet production process in Eq. 16, it is apparent that contractions of the virtual photon effective vertex given by Eq. 17 with δ µ− will give the result to be zero. This is a consequence of our choice of gauge. We therefore do not have any UV divergences in the (V 1) amplitude. There are however singularities arising from the l − 3 = 0 pole in the gluon propagator. They will be regulated by imposing a lower cutoff in the integral over z l at z 0 = Λ − 0 /q − , where Λ − 0 was specified in Eq. 4. To compute the terms proportional to logarithms in z 0 we will extract the contributions in Eq. 145 that are proportional to 1/z l . In terms of the constituent integrals given in Appendix B that appear in the computation of virtual graphs, where R (V 1) To obtain this expression, we redefined l 1⊥ − l 2⊥ → l 1⊥ in Eq. 145. The constituent integrals here have the (finite) expressions in d = 2 dimensions, where α 1 and α 2 are Feynman parameters. For the process (V 1) the arguments V ⊥ and ∆ are defined in terms of the factors (after imposing l 1⊥ − l 2⊥ → l 1⊥ ) in Eq. 148 as Interestingly, these arguments V ⊥ and ∆ appearing in the constituent integrals too can be expanded in terms of the gluon loop momentum fraction z l as where the coefficients c j (j = 1, . . . , 5 ) are different for each process. Using the above forms of the arguments one can perform the integration over z l in Eq. 150 to extract the rapidity divergent term and a remainder piece which is finite.
In general, we can write the amplitude for each processes as where The first piece is obtained from the contributions proportional to 1/z l . For the process (V 1), this contribution is the one given by Eq. 151. Following the same logic as articulated for the self-energy computations, it can be written (after extracting the rapidity logarithms) as the sum of a divergent contribution and a remainder term: R µα . The remainder term is comprised of terms in the amplitude which are proportional to 1/z l but are devoid of logarithms in z 0 . We can combine this piece with other finite terms in the amplitude that are not proportional to 1/z l . The contribution of the latter is denoted above by R (V β) (II);µα and is computed in Appendix I 4 for the different processes.
With this decomposition, we can finally write the net amplitude from the processes in Fig. 26 as In particular, for (V 1) we have where the divergent term is The integrals I v;log appearing in the above equation are the components of the constituent integrals in Eq. 152 that are proportional to logarithms in z 0 . These depend only on the coefficients c 1⊥ and c 3 of the arguments V ⊥ and ∆ of the constituent integrals when they are decomposed in the form shown in Eq. 154. This is a very general feature of our computation and we will express the divergent pieces for the other processes in terms of these integrals. Their expressions can be obtained as where repeated indices are summed over. For (V 1), the coefficients c 1 and c 3 appearing in these integrals are respectively Using the same techniques as described previously, it is straightforward to show that the UV divergent parts of the amplitude for the remaining processes (V 2), . . . , (V 6) are proportional to δ µ− and hence yield zero. In addition, there are no rapidity divergent contributions from the processes (V 2) and (V 5). Therefore we need to compute only the finite pieces R (V β) (II);µα (β = 2, 5) for these diagrams. The divergent term for (V 6) is similar to the one for (V 1). The computation of the divergent terms for (V 3) and (V 4) is especially tedious but can be carried out along the same lines as discussed here. These computations are cumbersome because we encounter more denominators whose evaluation, in turn, requires additional Feynman parameters in the transverse momentum integrations. However the arguments in these constituent integrals can again be written generically in terms of the gluon loop momentum fraction z l , as in Eq. 154. We will provide the results for the divergent pieces of (V 4) in Appendix G. The expressions for (V 3) can be obtained following similar methods.
Finally, the finite contributions from these 6 graphs ((V 1) − (V 6)) can be written compactly as The detailed computation of each of these terms is provided in Appendix I 4.

Contributions to T
(2) V : These are the processes labeled (V 7) − (V 12) which are obtained by interchanging the quark and antiquark lines in the diagrams shown in Fig. 26. As in Eq. 139 their net amplitude can be written as where In order to obtain T V , we have to impose the replacements defined in Eq. 66 along with a change of sign in the various terms constituting T (1) V in Eq. 140.

D. Vertex graphs with free gluon propagator
There are two broad categories of processes with six diagrams each that have the topology of vertex corrections with the gluon not crossing the shockwave. Three representative diagrams belonging to each of these classes are shown respectively in Figs. 28 and 31. The other half can be obtained by interchanging the quark and antiquark lines. The second class of processes labeled (V 19) − (V 24) actually resemble final state interactions between the quark and antiquark but we will also include them in this section. In the following sections, we will demonstrate the computation of such amplitudes by considering one representative diagram from each class, with the full result given in the Appendices.

Contributions to T
is the sum of the perturbative contributions from each process which can be decomposed into a divergent and a finite piece. We show below the computation of the divergent pieces constituting by considering the representative diagrams (V 13) and (V 15). The latter has a nested photon in the loop and hence has a topological similarity to the graph (S24) in Fig. 24. We will see that there are indeed divergent pieces that arise from this graph which cancel the divergences appearing in the third line of Eq. 138.
The amplitude for this process is given by The function Γ µ (l + 1 , l − 1 ) represents the gluon loop contribution to the γ * qq vertex. Unlike the processes discussed in the previous sections, where the sign of the loop longitudinal momentum l − 2 is determined by the location of the shock wave, there is nothing preventing l − 2 from taking both positive and negative signs. The expressions for the loop contribution for these two cases are derived in Appendix F and are given respectively by Eqs. F3 and F5. Using the delta functions present in the expressions for the effective quark vertices in Eq. 26, it is trivial to perform the l − 1 integration. We get l − 1 = −p − and the overall momentum conserving delta function . In order to perform the l + 1 integration, we need to separately analyze the two cases discussed above because the location of the l + 1 poles are dependent on the sign of l − 2 . The contribution to the amplitude for the process will then be the sum of the individual results obtained for these two cases.
• Case A: For 0 < l − 2 < k − + k − γ we have the following locations of the l + 1 poles, The second pole comes from one of the denominators in Eq. F3. We will enclose the pole at l + 1 | a for convenience.
• Case B: For 0 > l − 2 > −p − we have the following locations of the l + 1 poles, In this case, we will enclose the pole at l + 1 | c .
Using these results in Eq. 170 and finally subtracting the "no scattering" contribution we can write the amplitude for (V 13) as where and In this expression, we defined z l = l − 2 /q − and introduced and If we use the identity in Eq. D3 to expand the terms in the r.h.s of Eqs. 174 and 175 we can infer that the singular logarithms in rapidity will come from terms proportional to 1/z l whereas the UV divergent pieces will come from terms proportional to l 2 2⊥ l i , l i 2 l j 2 and l 2 2⊥ . We regulate these divergences using dimensional regularization in d = 2 − dimensions. We will therefore require the following constituent integrals to compute the divergent pieces of the vertex correction terms in Eqs. 174 and 175, where and α is a Feynman parameter. The forms of these arguments are different for the two cases. Employing all the above results, and a fair amount of algebra involving gamma matrices, we finally arrive at the following result for the amplitude (V 13): where In an identical manner, the amplitude for (V 14) is computed to be finite;µα (l 1⊥ ) . (182) One can therefore rewrite the amplitudes for (V 13) and (V 14) in terms of the amplitude for the leading order processes (labeled LO : (1, 2) in Fig. 8) as, We will now outline the computation of the amplitude of the graph (V 15). This is considerably more tedious than the previous two diagrams because the presence of the real photon in the loop allows for two possible regimes for the integration over l − 2 when l − 2 > 0, namely To see this explicitly, we will start with the amplitude for (V 15) (see Fig. 30 for the Feynman diagram with momentum assignments) given by As usual, the integration over l − 1 can be trivially performed using the Dirac delta functions contained in the effective vertices. The integration over l + 1 and l + 2 will be performed using the theorem of residues where the choice of contours is decided by the sign and magnitude of l − 2 . Expanding the denominators of the propagators we can see that there are three l + 2 poles located at Clearly the locations of l + 2 | a and l + 2 | b are well defined by the signs of the external longitudinal momenta but the pole at l + 2 | c can be below or above the real l + 2 axis respectively depending on l − 2 being positive or negative. We therefore have two cases similar to the computations of the graphs (V 13) and (V 14) in Fig. 28. However for the case l 2 > 0 we can have two separate contributions depending on the magnitude of l − 2 relative to k − . For l − 2 > 0, if we enclose the pole at l + 2 | b in the contour integration over l + 2 we have four l + 1 poles located at For 0 < l − 2 < k − we have the poles at l + 1 | a and l + 1 | b located below the real l + 1 axis whereas the remaining two are located above. We will deform the contour clockwise in our computation. For k − < l − 2 < k − + k − γ we have the poles at l + 1 | a,c,d located above the real axis and that at l + 1 | b located below which we will enclose through a clockwise deformation of the contour. The total contribution to the amplitude for the case l − 2 > 0 is therefore obtained by summing these individual contributions. The case l − 2 < 0 is simpler because we have a single pole located above the real l + 1 axis whose contribution obtained through an anticlockwise deformation can be computed using Cauchy's residue theorem.
The UV divergent terms are obtained from pieces proportional to the constituent integrals I There are no rapidity divergent pieces from this graph. We can finally write the amplitude for (V 15) as where contains the divergent piece that exactly cancels the residual divergence in the sum of the contributions from (S23) and (S24) in Fig. 24.
The finite contributions to the amplitudes for these processes are part of where the expressions for R (V β) finite (β = 13, 14, 15) are provided in Appendix I 6. As in the previous cases considered, Eq. 66 allows us to obtain the corresponding expressions for the finite contributions of the q ↔q interchanged processes (V 16) − (V 18). The amplitude for these processes is given by

Contributions to T
We have again expressed the perturbative contributions from each process as the sum of a divergent and a finite part. We will see that there are only rapidity singularities associated with these processes. We will detail below the computation for one such representative diagram (V 19). The amplitude contributions for the other two diagrams in Fig. 31 can be obtained following similar methods. In order to obtain the expressions for the q ↔q interchanged counterparts of diagrams (V 19) − (V 21) we will use Eq. 66 in  Fig. 31 with momenta and directions shown. This is representative of final state interactions between the quark and antiquark after they scatter off the shock wave.
for this process is given by where the free quark and gluon propagators are given by Eqs. 53 and the effective vertex for the dressed quark propagator is given by Eq. 26.
The integration over l − 1 can be trivially performed using the delta functions contained in the effective vertex factors resulting in an overall longitudinal momentum conserving delta function . We next perform the contour integration over l + 1 by enclosing the single pole situated below the real l + 1 axis. This results in the following expression for the amplitude: where the numerator and denominator are respectively given by and .

(195)
If we use the identity in Eq. D3 for the numerator, we see that it is at most proportional to l + 2 and hence the contour integration is well defined. There are three l + 2 poles located at .
The location of the first pole clearly depends on the sign of l − 2 . We therefore have the following cases: • Case A: For 0 < l − 2 < k − + k − γ we have l + 2 | a and l + 2 | b located below the real l + 2 axis whereas l + 2 | c is above. We will deform the contour anticlockwise to enclose this pole.
• Case B: For 0 > l − 2 > −p − we will deform the contour clockwise to enclose the pole at l + 2 | b . We are only interested in extracting the divergent pieces for these two cases. Redefining l 1⊥ + l 2⊥ → l 1⊥ we can get rid of the transverse phase containing l 2⊥ in Eq. 193. The denominator is now proportional to l 6 2⊥ in the limit of large l 2⊥ . Further, using the fact that terms proportional to δ µ− contributes zero to the cross-section, we can easily see that (V 19) is UV finite. There are however rapidity singularities and to extract them we need to isolate the terms proportional to 1/z l . Following the conventions used in Sec. V C we can write the amplitude for (V 19) as where finite;µα (l 1⊥ ) .
We have here terms proportional to 1/z l for the two cases A and B denoted respectively by R . From the first set of terms, we will get a logarithmically divergent contribution leaving behind a finite remainder. Adding these remainders with the other finite contributions R (V 19);A,B (II) results in the net finite contribution from this process: where the remainders are given by With these definitions in mind, we will now write down the expressions for the contributions proportional to 1/z l for the two cases considered above.
where the coefficients multiplying the constituent integrals are given by The expressions for the constituent integrals are given in Eq. 152. We can always express the arguments V ⊥ and ∆ of these integrals in terms of the gluon loop momentum fraction z l , as in Eq. 154. For the purposes of the discussion in this section where we extract the logarithmic singularity from Eq. 201 we will only require the expressions for the coefficients c and c . We can now finally write the rapidity divergent piece for Case A as where all but one of the integrals 8 appearing above are given in Eqs. 161. The arguments appearing in these can be obtained for case A as where the coefficients R (1,2,4) are the same as in Eq. 202 and The divergent piece for Case B can now be obtained as with c 1⊥ and c 3 defined in Eq. 204.
Combining the results for the two cases, we obtain the net divergent contribution to (V 19) as are given respectively by Eqs. 203 and 207. The expressions for the divergent pieces of the other contributions are provided in Appendix H. 8 The only missing expression is for the integral I

VI. CONSTRUCTING THE INCLUSIVE PHOTON+2 JET CROSS-SECTION
In this section, we will define jets 9 using the small cone algorithm of Ivanov and Papa [82] to extract the collinear and soft contributions from the expressions for the squared amplitude in the relevant diagrams of real gluon emission processes. These jet definitions can then be employed for amplitudes with virtual loops. When their cross-sections are combined with gluon emission contributions including collinear and soft divergences, we will obtain a finite crosssection for inclusive photon production in association with a quark-antiquark dijet.
Following the general definitions in Sec. II, the differential cross-section for inclusive γ + qq production at NLO can be written as whereX NLO;parton µν represents the NLO contributions to the hadron tensor, in analogy to its LO counterpart in Eq. 34. Its computation was outlined in Eq. 42 and carried out in the previous two sections. The superscript "parton" indicates that the various components that build up the hadron tensor in Eq. 209 are all calculated at the parton level. We will now discuss how to promote these quantities to the level of jets and shall construct the inclusive cross-section for γ + 2 jets.
As a first step towards writing down our final result, we refer the reader back to Table II. This table contains all the elements to construct the cross-section organized in terms of their color structures. Those with identical structures are placed in the same row, which makes the cancellations of divergences transparent. We begin with the virtual graphs discussed at length in the previous section; the structure of their divergences is apparent at the amplitude level. For most of these graphs, the divergent pieces are proportional to the LO amplitude; hence their interference contributions with the LO amplitudes are proportional to the LO cross-section.
We observed that the UV divergences arising from the self-energy graphs with dressed gluon propagators cancel with the divergent contributions from the self-energy graphs (S13) − (S16) where the gluon does not scatter off the shock wave. A similar cancellation of divergences that are proportional to the LO cross-section also occurs between the graphs (S21) and (S22). In addition, there are a few virtual gluon exchange processes, unique to photon+dijet production, which do not have a divergence structure proportional to the LO cross-section. In these graphs, the photon is nested in the gluon loop for both self-energy and vertex contributions. These correspond to the processes (S24) in Fig. 24 and (V 15) in Fig. 28 (and their q ↔q counterparts). However as shown in section V B, there is an intricate cancellation of divergences that takes place between the net contribution of self-energy graphs (S23) and (S24) in Fig. 24 and the vertex contribution from (V 15) in Fig. 28. The rapidity divergences also cancel between the diagrams {(S21), (S22)} and {(S23), (S24)}; they therefore do not contribute to the JIMWLK kernel. From the net UV contributions of virtual graphs, we are therefore only left with the divergences from vertex contributions with free gluon namely (V 13) and (V 14) (and their q ↔q counterparts).
We can now add the amplitudes for all the allowed virtual graphs and obtain the following result for their contribution to the NLO hadron tensor: The parton level LO hadron tensor is defined in Eq. 35. Recall that · · · corresponds to the CGC averaging in Eq. 23. The operator H V here contains bilinear functional derivatives in the classical gauge field A + cl ; as we will discuss in the next section, the action of these on the color structure Ξ of the leading order cross-section given in section II generates the leading log color structures of the NLO cross-section. The computation of the various pieces that constitute the finite contribution from the virtual graphs (denoted above byX NLO;parton µν;finite virtual ) are spelled out over seven subsections in Appendix I.
In this section, we will focus on the divergent structures shown in the first line on the r.h.s. These are the residual collinear singularities that survive after cancellations of UV divergences between different individual graphs in the virtual amplitude. One must therefore combine the cross-section containing these divergences with those of the real qqγ + g cross-section to obtain further cancellations.
For real gluon emission processes, it is difficult to extract their soft and collinear structures at the amplitude level. One needs to evaluate the squared amplitude and then integrate over the phase space of the outgoing gluon to recover these. In the notation of Table II, contributions from processes in which the gluon interacts with the shock wave (denoted as T (1) R ) do not contain soft or collinear singularities. This is because the gluon gains a net transverse momentum from the cumulative effect of successive "kicks" received during multiple scattering off the shock wave. Likewise, any soft gluon emitted before the shock wave is not energetic enough to cross it without being reabsorbed. This is true for the squared amplitude proportional to T R as well as the interference contributions of such graphs with those in which the gluon does not scatter off the shock wave; the latter are proportional to T Table II. The only kind of divergences present in these contributions are the small x logarithmic singularities which can be isolated by taking the slow gluon limit, k − g → 0 in the results obtained in Sec. IV. We will demonstrate these in detail in the next section. One can then obtain the finite contribution to the cross-section from T where the term proportional to the small x logs can be expressed as the convolution of an operator H (1) R (analogous to H V in Eq. 210 above) acting on the leading order hadron tensor. This too will be discussed in the next section. The finite pieces can be evaluated numerically; doing so will be a topic for future work.
The gluon emission processes of interest in this section will be the ones in which the gluon does not cross the shock wave and we have a region of phase space where it can be soft or collinear with respect to the quark or antiquark. Depending on whether the gluon is emitted from the quark or antiquark, we have denoted the contributions from these kinds of processes respectively by T  Table II that we get different results for the respective color structures in the NLO cross-section depending on whether the processes constituting T (2,3) R interfere between themselves or with their q ↔q interchanged counterparts. This dependence on color flow at the loop level is absent for the case discussed in [54][55][56] where the qq is projected on to a singlet final state. We will show that these differing color structures have interesting implications for soft gluon factorization.
Since we are not integrating over the phase space of the quark and antiquark, combining the NLO cross-sections for qqγ and qqγ +g will not in general cure all infrared (IR) singularities. We will still have a remnant collinear divergence which can be absorbed in the jet fragmentation function of the quark or antiquark and this can be interpreted as contributing to the evolution in energy of the fragmentation function. Conversely, we can construct an IR-safe cross-section for photon+dijet production by using a jet algorithm that restricts the phase space for the final state gluon.
Following [55,82], we will work in the small cone approximation (SCA) in which the extent of the jets in the rapidity-azimuthal angle (Y, φ) plane is small, or more quantitatively, the jet cone radius R is not too large (R 2 1). In this approximation, one then systematically expands the partonic cross-sections around R = 0. The dependence on R is of the form A ln(R) + B + O(R 2 ), where the coefficients A and B can be evaluated analytically and the terms that are power suppressed in R 2 are neglected.
For the inclusive computation of the photon+2 jet cross-section there are three possible cases that we need to consider: • The gluon is inside the cone of either the quark or antiquark, • The gluon is outside the cone, • The gluon forms one of the jets, while the other jet is formed by a qq pair.
We will not consider the third sort of contribution in this list because it does not have a collinear divergence and is hence sub-dominant in the SCA. A short proof is provided in Appendix J.
We will therefore first isolate the singularities from the region of phase space where the real gluon is collinear to the quark or antiquark and shall use the SCA framework to identify when two partons form a jet. For a jet cone radius R (R 2 1), two partons i and k with respective momenta p i and p k will form a jet 'J' carrying a momentum equal to the sum of their momenta if both partons satisfy the condition Here ∆φ i,k is the difference of the azimuthal angles between parton i (k) and the jet; ∆Y i,k is the corresponding rapidity difference. In the SCA, the jet constituted of the partons i and k is considered on-shell (upto O(R) corrections) and hence its momentum can be written as The quantities on the l.h.s of the inequality in Eq. 212 are given by [82] Here z i,k = p − i,k /q − . We will introduce the "collinearity" variable which approaches zero when the partons i and k are collinear (p k⊥ → z k /z i p i⊥ ). It is then possible to rewrite and express the quantities in Eq. 214 in terms of C ik,⊥ . The small cone condition in Eq. 212 can then be equivalently written in terms of this collinearity variable as We will now use the above definitions to extract the collinearly divergent contributions from the processes in which a quark and gluon constitute the first jet 'J' and the antiquark constitutes the second jet 'K'. The parton level Feynman graphs for these contributions are shown in Fig. 20. The corresponding expression for the quark-antiquark interchanged diagrams are obtained simply by J ↔ K symmetrical replacements. If we look carefully at the diagrams in Fig. 20, processes labeled (R12) − (R15) have a gluon emitted by an on-shell quark with the remaining structure identical to those of leading order (LO) processes depicted in Fig. 8.
It is physically intuitive that the divergent structure of the amplitude squared for these processes will be proportional to the LO amplitude squared. We will see that this is indeed the case. For the diagram labeled (R11) (see Fig. 21) the collinearity of the gluon with respect to the quark is lost because the photon emitted after the gluon imparts a virtuality 2k.k γ to the quark. However in the soft photon limit k γ → 0, the quark goes on-shell and can be collinear to the gluon. In the case of the amplitudes for the graphs (R12) − (R15) in which the gluon is emitted after the photon, we see from Eqs. C10-C13 that there is a possible divergence in the collinear limit k g⊥ → z g /z q k ⊥ . Because there is a similar term in the numerator of these expressions, the structure of the singularity becomes transparent only at the level of the squared amplitude.
It is also easy to check using the expressions for the amplitude that for interference contributions between diagrams in which the gluon is emitted from the quark in the amplitude and from the antiquark in the conjugate amplitude (or vice versa) there are no divergent terms in the collinear limit. Collinear divergences arise only when the gluon is emitted and reabsorbed by the parton to which it is collinear.
For the case where quark and gluon form the jet J, and the antiquark constitutes jet K, we will introduce the jet variables, so that the d-dimensional phase space differential measure transforms as gluon pols.
for the sum over gluon polarizations, it is a straightforward exercise to show that the following relation holds in the collinear limit C qg,⊥ → 0, When we integrate over the phase space of the gluon, the small cone condition given by Eq. 217 restricts the integration over C qg,⊥ to be Here d = 2 − as previously. We will denote this collinearly divergent (dominant) contribution to the cross-section from the amplitude squared of processes (R12) − (R15) asX NLO µν;collinear . From Eq. 221, and employing the phase space factors on the r.h.s of Eq. 219, this is given bỹ whereX LO;jet µν is shorthand for the LO hadron tensor for the production of a photon plus the quark-antiquark jets J and K.
This quantity is equivalent to the LO hadron tensor for γ + qq production if we make the following replacements in Eq. 35: The lower cutoff on the gluon momentum fraction z g is set to z 0 = Λ − 0 /q − as in the previous sections. We will use the following result in dimensional regularization for the transverse integral, The expression for the collinear contribution to the cross-section from the q ↔q interchanged counterparts of (R12) − (R15) is obtained simply by J ↔ K symmetry. The net contribution from these processes is We observe clearly in the above expression that the collinear divergences cancel and so do the terms proportional to ln 2 (1/z 0 ). We are however left with a term from Eq. 226 that contains a soft-collinear divergence proportional to ln(1/z 0 ) ln(R) (in addition to finite pieces for a fixed cone size R).
Recall however that this result is obtained in the collinear limit C qg,⊥ → 0 of the real graphs as shown in Eq. 222. The dominant contributions to the inclusive cross-section in this limit are when the gluon is inside the cone of the quark or antiquark jet. When one additionally requires the gluon to be soft, there can be contributions where the gluon is outside the jet cone. As we will now demonstrate, the jet algorithm when applied to the case of a jet formed by a single parton generates soft-collinear terms with opposite sign that exactly cancel those obtained in the collinear limit.
In the soft gluon limit k g → 0, the contributions to the amplitude from the five processes in Fig. 20 and their q ↔q counterparts (given respectively by Eqs. 67 The sub-leading pieces come from the process (R11) (and its q ↔q counterpart) where the gluon is emitted from an internal fermion line as well as from the next-to-leading soft terms in the amplitudes for (R12) − (R15) (and their q ↔q counterparts) in the expansion around k g = 0.
Squaring this expression, summing over the gluon polarizations, and taking the CGC average over all possible static color charge configurations, the dominant soft gluon contribution to the NLO hadron tensor (and hence the NLO cross-section) isX where we have broken up the soft gluon contribution into two parts. The first contribution comes from the terms T Table II and possess a color structure similar to that of the LO hadron tensor in Eq. 35: This is however not true for the interference contributions T R +c.c in Table II; as noted previously, such interference contributions do not contain collinear divergences. These however have the following dominant contribution to the cross-section in the soft gluon limit 10 , where C LO µν is the LO coefficient function, There are no soft-collinear divergences in Eq. 231. The only divergent contributions are those in rapidity that are recovered when taking the slow gluon limit. We will therefore focus on extracting the soft gluon contributions to the jet cross-section from Eq. 230. Since the collinear contribution obtained in Eq. 226 already includes the soft-collinear limit, adding this contribution to the expression in Eq. 230 will result in a double counting of such divergent pieces. To avoid this, we will restrict the region of integration in the first term on the r.h.s of Eq. 230 to ensure the emitted gluon is outside the cone of the jets formed by the quark/antiquark: where C 2 qg,⊥ | min. = R 2 p 2 J⊥ /z 2 J for the first integral is obtained by imposing the small cone condition in Eq. 212 on a jet formed by a single parton [82]. The upper bound for the integral over the collinearity variable can in principle be infinity. Hence the integration can be performed using dimensional regularization and we will obtain a similar result as in Eq. 225. We get finally, As promised, we see that the terms in the first line on the r.h.s are identical but have the opposite sign to the corresponding terms in Eq. 227, therefore canceling the final remaining collinear divergence. The double log appearing in the second line of the equation above is contained in the "slow" gluon limit (k − g → 0) of our results. As we noted previously from our discussion of similar divergent terms in the virtual graphs, this is the double log limit of the BFKL equation. We are therefore double counting here because the soft gluon sector is a subset of the slow limit. In order to obtain finite contributions from the squared real amplitude, we must subtract the pieces from the soft gluon limit that contribute to small x evolution 11 .
More specifically we will absorb the kernels obtained in the slow gluon limit (of the real unscattered gluon contributions constituting T R → H R , whose structure we shall discuss further shortly in the next section. When we combine this sum with the result from Eq. 227 we obtain finallỹ We see that the soft-collinear pieces cancel leaving terms that have the generic form A ln(R) + B as expected in the SCA. As discussed above the divergent pieces coming from the slow gluon limit that are contained in the soft gluon limit can be isolated and absorbed in a modification of H R to form the operator H R . As we will discuss further in Sec. VII, H V + H R = H LO , where the r.h.s is the LO JIMWLK Hamiltonian that we introduced previously in the introduction. Finally, the finite contributions appearing in the second line together with the first term on the r.h.s of Eq. 235 constitute the NLO impact factor for inclusive photon+dijet production. The computation of these terms is the principal result of this work.
We conclude this section with a few comments. Firstly, we observe that there are no Sudakov double logs in our computation-this also holds for inclusive dijet production case in the soft photon limit of our result. These logs appear due to the lack of complete real-virtual cancellations arising from a constraint imposed on the process, an example being dijet production in back-to-back correlations 12 . Since we do not impose any such kinematic constraints, this explains our observation.
Secondly, soft gluon factorization is violated. This factorization, in the soft gluon limit of k g → 0, corresponds to a convolution of the LO dijet cross-section with the well known eikonal kinematic factor |k β /(k.k g ) − p β /(p.k g )| 2 , where the latter is contracted with the sum over gluon polarizations. This violation is a consequence of the differing topologies of the color structures that contribute towards soft and collinear divergences. However we can rewrite the second term on the r.h.s of Eq. 229 (which violates the soft gluon theorem) as a term that obeys soft gluon factorization and another that violates factorization. With this, we can rewrite Eq. 229 as The dipole (D) and quadrupole (Q) operators defined in Eq. 39 correspond to particular boundary conditions at x − = ±∞. There have been recent developments that relate soft gluon theorems to the existence of infinite dimensional so-called BMS symmetries [98] and to a color memory effect in Yang-Mills theory [99]. A dictionary between this language in the Regge limit of QCD and that of the CGC was established in [100], and involves identifying the spacetime rapidity η = ln(x − ) in the latter with retarded time in the former. An interesting question is whether the soft gluon theorem can be restored by requiring that the factorization violating structure (Q − DD) vanishes by a modification of boundary conditions at x − = ±∞. This may be equivalent to defining asymptotic states/propagators that project out this color structure at x − = ±∞. We will return to this interesting topic in future.

VII. HIGH ENERGY LEADING LOG RESUMMATION
In this section, we will consider the "slow" (relative to the virtual photon LC momentum q − ) gluon limit of our results for real (z g → 0) and virtual (z l → 0) diagrams. This will allow us to isolate the soft divergences in rapidity; we will show explicitly that these terms provide a nontrivial derivation, in the evolution of the projectile, of the JIMWLK renormalization group equation.
We begin by examining the amplitudes for virtual processes in the z l → 0 limit, and subsequently, real emission graphs. For the self-energy contributions with the dressed gluon propagator, consider the process (S1) given by Eq. 88: where the parameters v (S1) 1,2⊥ and ∆ (S1) 1,2 are all proportional to z l as can be seen from Eq. 89. In the slow gluon limit z l → 0, this expression simplifies greatly to give where the contribution R LO:(1) µα (l 1⊥ ) to the perturbative amplitude at LO was defined previously in Eq. 29. The integrals over l 2,3⊥ represent the spatial derivative of the two dimensional free gluon propagator. We will use its coordinate space expression 13 where the r.h.s represents the well-known Weizsäcker-Williams field created by the boosted quark-antiquark pair. Substituting this in Eq. 238 then gives There is a logarithmic divergence in the kernel for x ⊥ → z ⊥ under which the color structure also reduces to that for the LO process times the Casimir C F . This is precisely the 1/ singularity multiplying ln(1/z 0 ) that appears in the momentum space result for (S1) given in Eq. 78. Note that the color structure there also has the form C F (Ũ (x ⊥ )Ũ † (y ⊥ ) − 1) for these divergent pieces. These divergences cancel in observables reflecting the fact that there are no divergences in the LO JIMWLK Hamiltonian in the limit x ⊥ , y ⊥ → z ⊥ . One can show in a similar fashion that (S2) and (S3) vanish in the z l → 0 limit. The three remaining processes have an identical structure to Eq. 240 albeit with different R-functions corresponding to the structure of their LO counterparts. Combining these results, the six contributions in Fig. 22 give The quark-antiquark interchanged diagrams are obtained by replacing x ⊥ ↔ y ⊥ in this equation. We next consider the slow gluon limit of the self-energy corrections containing a free gluon propagator. Half of these 24 processes are depicted in Fig. 24, with the other half obtained by interchanging the quark-antiquark lines. In the limit z l → 0, the loop contribution shown in Fig. 33) for the quark self-energy reduces from the expression in Eq. D7 to the simpler expression, Employing the identity we can deduce an identical form for this limiting expression as in the case of contributions with dressed gluon propagators. Using the above simplification, and the identity in Eq. 239, the amplitude for (S13) given in Eq. 111 simplifies to The above expression has the same form as Eq. 240 albeit with a negative sign and a different color structure. The corresponding expression for its quark-antiquark interchanged counterpart is obtained by imposing x ⊥ ↔ y ⊥ . A similar simplification occurs for (S14)-(S16) with the replacement R LO:(1) µα with the corresponding LO R-functions for these graphs.
Recall that as discussed previously in Section V B, the diagrams labeled (S17)-(S20) in Fig. 24 vanish in general. More nontrivial is that fact that the limiting forms of the last four processes labeled (S21) − (S24) in Fig. 24 (and their q ↔q interchanged counterparts) cancel amongst each another. This can be shown explicitly by using Eq. 242 for (S21) and (S23) in Fig. 24 and the simplification ofΣ α (k f , k γ ) in Eq. E4, for the graphs (S22) and (S24) shown in Fig. 34. Thus (S13)-(S16) (and their quark-antiquark exchanged counterpart) are the only surviving self-energy contributions with the free gluon propagator; the net result in the slow gluon limit of the 24 self-energy diagrams is simply We used here the relation (l ⊥ ) given in Eq. 33. We now turn to the structure of vertex corrections in the slow gluon limit. As an example, consider the amplitude (V 1) we worked out previously for the dressed gluon propagator, the results for which are spelled out in Eqs. 145, 146 and 147. After further redefining l 1⊥ − l 2⊥ → l 1⊥ , and using Eq. 239, we obtain in the z l → 0 limit, Similarly, it can be shown that the limiting forms for processes (V 2) and (V 5) are zero. Combining the contributions for the remaining three processes, we can write the total amplitude for the 6 processes in Fig. 26 under the z l → 0 limit as lim z l →0 The structure of this kernel is different from that of the self-energy corrections because here the Weizsäcker-Williams fields are emitted by both the quark and the antiquark. Since the kernel is symmetric under x ⊥ ↔ y ⊥ , the limiting expression for the amplitude under quark-antiquark interchange has the same form; note however that for the latter, the color structure is instead ( Likewise, one can show that the corresponding virtual amplitudes with free gluon propagators (in Figs. 28 and 31) are At the amplitude squared level, we have to consider the interference contribution of these virtual graphs with the leading order amplitude result in Eq. 32. As outlined in the introduction, we then need to perform a CGC averaging over all possible source charge configurations ρ A . Using the expressions for the loop contributions derived above in III. LLx structure of the amplitude squared from interference contributions of virtual graphs with LO processes. The kernels for the complex conjugates of the squared amplitudes in rows 3 and 4 are obtained by replacing x ⊥ → x ⊥ , y ⊥ → y ⊥ . The color structure and kernels for the complex conjugates of rows 1 and 2 are obtained by permutations of the coordinates.
the z l → 0 limit, and the color structures given Table II, we can obtain the leading logarithmic singular structure (in rapidity) of the squared amplitude for virtual graphs. These are summarized in Table III.
The coefficient function C LO in this table was defined previously in Eq. 232. The "CGC averaged" leading order squared amplitude constituting the LO hadron tensor in Eq. 22 can be expressed as with Ξ(x ⊥ , y ⊥ ; y ⊥ , x ⊥ ) = 1−D xy −D y x +Q xy;y x , where the dipole (D) and quadrupole (Q) traces over the lightlike Wilson lines were defined previously in Eq. 39.
We will now repeat the above exercise by taking the slow gluon limit of gluon emission diagrams. We will first take the z g (= k − g /q − ) → 0 limit in these amplitudes, take their modulus squared, and then integrate over the phase space of the emitted gluon. We begin by considering the contribution from the 10 graphs in which the gluon crosses the nuclear shock wave. (One half of these are shown in Fig. 18 and the other half are obtained via quark-antiquark exchange.) Recall that their combined contribution was expressed in Eq. 49 as where z r tot = z q + zq + z γ + z g and T R is defined to be As an example, consider the process (R1) for which the R-function (in Eq. 62) is where the I functions are proportional to modified Bessel functions of the second kind and are given in Eq. 63. The factors in the arguments of these functions multiplying which vanish in the z g → 0 limit. Under these conditions, we have and One can similarly derive limiting expressions for the remaining 4 diagrams in Fig. 18. The coordinate space structure of the kernel remains the same for these diagrams but they are proportional to different R LO 's. For the quark↔ antiquark interchanged diagrams, we can simply replace x ⊥ ↔ y ⊥ along with an overall change of sign to get the corresponding expressions in the z g → 0 limit. We obtain the z g → 0 limit of the amplitude from these 10 contributions as As a final demonstration, we will consider (R11), a diagram where the gluon is emitted by the quark and does not scatter off the background classical field. Its amplitude was given previously in Eq. 70. In the z g → 0 limit, this amplitude reduces to The limiting expression for the net contribution from the 5 processes in Fig. 20 therefore can be written as Similarly, the contribution from their quark↔antiquark interchanged counterparts are, These three amplitude structures therefore give 9 contributions at the level of the cross-section, many of which are Hermitian conjugates of each other. Since the final state process of interest is inclusive photon+dijet production, we will integrate over the phase space of the emitted gluon. To obtain coordinate space kernels, we will liberally use Eq. 239 and the relation between the two dimensional free propagator in momentum and coordinate space: We can now extract the leading logarithm in x (LLx) structures of the various real emission contributions to the differential cross-section at NLO. These are summarized in Table IV. The quantities in the first column represent the CGC averaged squared amplitudes integrated over the phase space of the emitted gluon, We now have all the essential elements we promised in the introduction of the paper to derive the JIMWLK evolution equation. Using the leading logarithmic structures summarized in Tables III and IV for virtual and real emissions we can organize the CGC averaged squared amplitudes (and their Hermitian conjugates) in a basis spanned by dipole, D and quadrupole, Q Wilson line correlators and their products DD and DQ. (For a general introduction to such basis structures, and how to compute them, we refer the reader to [101].) After some algebra, and use of the identity we can derive the following leading logarithmic structure for the CGC averaged amplitude squared at NLO In the above equation, z f is the scale of the longitudinal momentum of the virtual photon probe. It is now a straightforward but tedious exercise to show that the set of terms appearing within curly brackets in the above equation can be generated by the action of the LO JIMWLK Hamiltonian on the leading order color structure Ξ(x ⊥ , y ⊥ ; y ⊥ , x ⊥ ). The former is defined as [10][11][12][13][14] where The Wilson lines appearing above are in the adjoint representation of SU (N c ) and written in terms of A +,a (x − , x ⊥ ) as where we have defined A +,a (x ⊥ ) = +∞ −∞ dz − A +,a cl (z − , x ⊥ ). As first outlined in the JIMWLK papers, the proof involves extensive use of the identities, and Eq. 48 relating the fundamental and adjoint Wilson lines. Now using Eq. 251 for the CGC averaged amplitude squared at LO, one can write Eq. 265 as While we derived the l.h.s of this identity explicitly here for inclusive photon+dijet production, all the necessary elements to derive the r.h.s were obtained previously in [102].
In the CGC EFT, the expectation value of an operatorÔ is defined to be whereÔ[ρ A ] is the expression for the operator for a given charge configuration ρ A and W z [ρ A ] is a stochastic weight functional describing the probability density of that charge configuration at a momentum (fraction) z. If we now consider the result for the CGC averaged squared amplitude (or differential cross-section) upto NLO+LLx accuracy we can easily see that the leading logarithmic pieces can be absorbed in the weight functional to redefine the EFT at the evolved scale z f as In this equation, the functional dependence on ρ A enters through the dipole and quadrupole Wilson line correlators contained in Ξ(x ⊥ , y ⊥ ; y ⊥ , x ⊥ ) appearing in the LO amplitude squared.
Since the lepton tensor, and other prefactors remain the same, we can extend these arguments to obtain the JIMWLK evolution of the triple differential cross-section for photon+dijet production in small x DIS: Our proof is in the spirit of the JIMWLK derivation from the projectile side in [103] but is obtained by computing the full cross-section and then taking the slow gluon limit.

VIII. SUMMARY AND OUTLOOK
We presented in this paper the first computation of the NLO impact factor for the inclusive photon+dijet production in high energy electron-nucleus collisions. The triple differential cross-section for this process can be expressed as where L µν is the lepton tensor defined in Eq. 21 and In this expression, the finite termsX NLO;jet µν;finite are of order α S relative to the leading term. The explicit results we derived for these are the principal results of this paper. In order to obtain their NLO expressions, we showed in sections V and VI that one has to first isolate the ultraviolet, collinear and soft divergences respectively in the real gluon emission and virtual self-energy+vertex correction diagrams. For the virtual graphs containing gluon loops, as discussed in section VI, several already cancel between different contributions to the amplitude. Of the rest, the collinear and soft-collinear divergences cancel between the real and virtual graphs at the level of the squared amplitude.
Our treatment of these in Section VI necessitated introducing a jet algorithm; the small cone approximation (SCA) framework, corresponding to a small jet cone radius R, is particularly convenient for our task. The term proportional to α S ln(R) in Eq. 275 is an O(1) remnant of this procedure. The soft gluon limit (k g → 0) is a subset of the slow gluon limit of k − g → 0 but finite k ⊥ . The large logarithms in x arise in the latter; the resummation of these O(α S ln(1/x)) terms is discussed at length in section VII. Our discussion there provides a nontrivial derivation of the JIMWLK equation. Though this equation has been derived in a variety of ways, it is interesting to see it arise from an explicit Feynman diagram computation of a nontrivial final state.
An apparently technical point which is however of general interest is our observation in section VI that the JIMWLK kernel already contains pieces of what we isolated as soft-collinear divergences. These have to be subtracted from the jet cross-section to avoid double counting when the NLO impact factor is combined with small x evolution. The fact that slow gluon emission outside the jet cone satisfies JIMWLK evolution is consistent with this being a feature of nonglobal logarithms in jet physics explored by Banfi, Marchesini and Smye [104], identified with BK/JIMWLK evolution by Marchesini and Mueller [105] as well as by Weigert [106], and subsequently significantly developed by Hatta and collaborators in a number of papers [107][108][109][110]. (See also [111] for a recent discussion of this correspondence.) Our NLO computation of photons+dijets in DIS therefore combines JIMWLK evolution in both the spacelike evolution of the DIS wavefunction and the timelike evolution of dijets in the final state. This spacelike-timelike correspondence in A − = 0 gauge was noted previously by Mueller [96] and and is a quantitative implementation of a proposal in his paper.
The JIMWLK evolution equation describes the small x evolution of the gauge invariant weight functional W LLx xBj [ρ A ] which resums leading logs O((α S ln(1/x)) n ) and power corrections 14 O((Q s /Q) n ). An important development is that the NLO JIMWLK evolution equation that resums terms of order O((α 2 S ln(1/x)) n ) is now available (in addition to a significant body of work on the NLO BK equation). If we take advantage of these developments, we can promote The finite termsX NLO;jet µν;finite in Eq. 275 from the virtual loop contributions are given in Appendix I. The finite terms from the real gluon emission contributions to the cross-section are obtained by taking the modulus squared of the amplitudes in Eqs. 49, 67 and 73, integrating over the gluon phase space with a cutoff, implementing the SCA in section VI, and subsequently subtracting the pieces that contribute to leading log JIMWLK evolution.
The eventual goal of this project is to provide quantitative predictions for a future Electron-Ion Collider (EIC). As noted in the introduction, the computation of the finite piecesX NLO;jet µν;finite along with NLO BK/JIMWLK evolution, provide the necessary ingredients to compute photon+dijet production (and the associated measurement channels we identified) in e+A DIS to O(α 3 S ln(1/x)) accuracy. For the x reach of an EIC, this corresponds to an accuracy of O(α 2 S ) or ∼ 10%. This level of accuracy in such differential measurements is likely sufficient for the unambiguous discovery of gluon saturation.
The realization of this numerical project while clearly feasible is nevertheless a formidable challenge on several fronts. Firstly, the number and complexity of the finite contributions to the NLO impact factor are far greater than comparable expressions in the collinear factorization framework. This is because, unlike the latter, the large k ⊥ ∼ Q S from the target flowing through the diagrams generalizes functional forms in the collinear framework to nontrivial integrals that in many instances have to be performed numerically.
Further, the results obtained for the NLO JIMWLK Hamiltonian are not yet ripe for numerical evaluation. In addition to the sheer complexity of the NLO evolution kernels, there are subtle conceptual issues, first identified by Salam [25], regarding the regularization of these kernels-recent work in this direction, and references to the extant literature, can be found in [27,112]. A self-consistent treatment of NLO JIMWLK in our framework can be obtained by computing the leading ln(1/x) contributions to the next-to-next-to-leading order coefficient function for dijet production in e+A DIS. While challenging, the developments in this paper suggest it can be achieved on the required time scales.
Both the numerical and analytical developments suggested here are however beyond the scope of this work and will be pursued in parallel in future. As a final note, the methods, computations, and principal results of this paper are summarized in an accompanying letter [88].

Appendix A: Notations and conventions
The metric used in this paper is the −2 metric,ĝ = diag(+1, −1, −1, −1), where the 'hat' denotes quantities in usual spacetime coordinates. The light cone coordinates are defined as with the transverse coordinates remaining the same. The same definition holds for the gamma matrices γ + and γ − with the Dirac algebra given by where g +− = g −+ = 1 and g ij = −δ ij (i, j = 1, 2) are the nonzero entries of the metric tensor. In this convention,

Appendix B: Constituent integrals in computations of gluon emission and virtual gluon loop diagrams
In this Appendix, we will derive a generic expression for the tensor integrals that arise in the amplitude (and squared amplitude) computation of real emission diagrams and use it to extract some specific results. By taking a specific limit of this result, we will also get the expression for the tensor integrals that arise in the amplitude computation of virtual gluon exchange diagrams.
• Generic integral for gluon emissions: The most general tensor integral in d dimensions appearing in real gluon emission computations has the form The denominators D 1 , · · · , D n appearing on the r.h.s of the above equation have the form D i = (l ⊥ +v i⊥ ) 2 +∆ i . By introducing n Feynman parameters α 1 , · · · , α n , we can reexpress these integrals as By integrating out one of the Feynman parameters, it is always possible to write the denominator in the above expression as (l ⊥ + V ⊥ ) 2 + ∆, where V ⊥ and ∆ can be written in terms of the various v i⊥ 's and ∆ i 's with the Feynman parameters acting as coefficients. This will become more clear when we present the expressions for some specific cases below. The integral in Eq. B1 can therefore be written as The remaining Feynman parameters satisfy the condition To reduce this above integral, we first make the observation that Performing this operation p number of times, we can derive the identity • Generic integral for virtual gluon exchange diagrams: The most general tensor integral appearing in the computation of virtual graphs has the form Clearly this can be obtained in the r ⊥ → 0 limit of the integral in Eq. B1. We have provided explicit expressions for such integrals wherever they appear in the main body of the paper. To derive such expressions from the result given by Eq. B8 one needs to carefully take the limit r ⊥ → 0 in the MacDonald functions K ν ( r 2 ⊥ ∆) that will appear in the computation. The limiting result is given by (see Eq. 10.30.2 of [113]) Appendix C: R-functions appearing in gluon emission amplitudes In this section, we will present expressions for the R-functions (not given in the main text) embedded in the final expression for the NLO amplitude for gluon emissions: Here we defined T (1) 33. Quark self-energy diagram for quark momentum k f . This block can appear either in the internal or external legs. where It is a trivial exercise to show that where (. . .) represents the terms in between the gamma matrices in Eq. D1. Using this identity, we encounter the integrals over l + , While the first integral can be done trivially using contour integration, and Cauchy's residue theorem, the second integral should be examined more carefully. This is because unlike the first case where the contribution on the semicircle of the contour vanishes as the radius of the semicircle approaches infinity, it does not for integrals of the second kind. Therefore to obtain the value on the real line (which is the second integral in Eq. D4), we must subtract the semicircular contribution from the value obtained using the residue theorem. The location of the l + poles given by Eq. D2 depend on the sign of l − . We need not worry about the l − = 0 case because that will be regulated by imposing a cutoff at the value Λ − 0 . For 0 < l − < k − f the poles are on opposite sides of the real l + axis whereas for l − < 0 they are on the same side. Let us discuss the latter case first. When both poles are above the real l + axis, integrals of the first kind in Eq. D4 give a null result because the contour can be closed in the other direction without enclosing any pole. For integrals of the second kind, if we choose to close the contour above such that both the poles are enclosed, then we have the following result after subtracting the contribution from the semicircle (as the radius approaches infinity).
We will obtain the same result if we close the contour clockwise in which case the sole contribution will be the one from the semicircle. We will come back to this result later. Now for the case 0 < l − < k − f if we choose to close the contour clockwise enclosing the pole at l + a we get the following result for the second integral, where the quantity of interest is The choice of arguments for Γ µ will become clear as we proceed with the calculation. We will first perform the integration over l + 2 using complex contour integration. Using the identity in Eq. D3 and the fact that terms proportional to γ − γ µ γ − = 2δ µ− γ − yield zero after contraction with the intermediate photon vertex in Eq. 17, we see that the numerator is at most proportional to l + 2 whereas the denominator is proportional to (l + 2 ) 3 . As such, the contour integration can be performed without any additional complications. The three l + 2 poles are, where l + 2 | b and l + 2 | c are located below and above the real l + 2 axis respectively for any l − 2 . The location of l + 2 | a however depends on the sign of l − 2 . In the following, we will obtain the generic expressions for Γ µ for two cases • Case A: For 0 < l − 2 < q − + l − 1 we have the poles at l + 2 | a and l + 2 | b located below the real l + 2 axis whereas l + 2 | c is above. We will therefore deform the contour anticlockwise to enclose the pole at l + 2 | c . Using Cauchy's residue theorem and making the momentum redefinition l 2⊥ − l 1⊥ → l 2⊥ with l 1⊥ remaining unchanged we obtain the following expression for Γ µ where The choice of l + 1 and l − 1 as the arguments of Γ µ is now transparent from the above expressions. From the identity derived in Eq. D3, it is clear that the numerator in Eq. F3 above has to be evaluated at l + 2 | c which will have a residual l + 1 dependence. The value of this will depend on the choice of contour taken for the l + 1 integration for the full computation. The same is also true for the factor ∆ V 1 given by Eq. F4.
• Case B: For 0 > l − 2 > l − 1 we have the poles at l + 2 | a and l + 2 | c above the real l + 2 axis while the pole at l + 2 | b is below. For convenience we therefore choose to deform the contour clockwise. Performing the same momentum redefinition l 2⊥ − l 1⊥ → l 2⊥ as above and applying Cauchy's residue theorem we get where Once again the above expressions are to be evaluated at the l + 1 pole enclosed by our contour of choice in the full computation. In this section, we will discuss the strategy to compute the remaining divergent contributions (those not provided in the main text) that constitute T (1) V in the contributions to the amplitude from the six diagrams in Fig. 26. Recall that this amplitude can be expressed as and each R (V β) , β = 1, . . . , 6, can be decomposed as the sum of a divergent and a finite part as In Sec. V C, we demonstrated that there are no UV divergent contributions from these diagrams for our choice of gauge. We are only left with singular pieces in rapidity for four out of the six allowed processes in this category. For (V 2) and (V 5), we have purely finite contributions; we can therefore write, The expressions for these finite pieces will be given in Appendix I 4.
We begin with the process (V 6) which has a similar structure to the process (V 1). We can write where the divergent contribution is obtained as The integrals appearing above are given by Eqs. 161. For (V 6) we obtain the following expressions for the arguments c 1⊥ and c 3 of these integrals, We will now present the results for (V 4) which has a similar structure to that of (V 3). We will not provide expressions for the divergent pieces of (V 3) here because they are similar to that of (V 4) but more lengthy. The latter is due to the fact that in the amplitude computation of (V 3) the contour for the integration over l + 1 encloses two poles on the same side of the real axis. There are therefore two separate contributions which need to be added in order to obtain the final divergent piece. The divergent terms computed for these processes contribute towards the leading logarithmic evolution of the LO result and can be absorbed in a redefinition of the weight functional W Λ − 0 [ρ A ] describing color sources as described in the introduction. This is also explicitly shown in the JIMWLK derivation discussed at length in section VII.
For (V 4), we can write where the divergent pieces for (V 4) are obtained from terms in the amplitude that are proportional to 1/z l . These specific terms can be written as The finite expressions for the constituent integrals appearing in the above equation are as follows: The remaining integrals can be obtained by equating two of the indices in the expressions provided above. The arguments V ⊥ and ∆ can always be expressed in terms of z l as shown in Eq. 154. The integration over z l can now be performed and we can express the divergent piece proportional to logs in z 0 in terms of integrals over the Feynman parameters. These integrals only depend on the coefficients c 1⊥ and c 3 appearing in Eq. 154 whose individual expressions in turn depend on the process of interest. The divergent term in (V 4) can now be finally written as ) + 4z q q − δ iµ γ α γ j γ − − 2q − γ i γ µ γ α γ j γ − − 2zqq − γ i γ µ γ j γ α γ − × I  ) .
The expressions for some of the integrals appearing in the above equation are provided below. The rest of them can be obtained by putting i = j in the expressions given.
The finite expressions for these constituent integrals are given in Eqs. G10. It should also be noted that the limits of integration over z l are different for the two cases we are considering. For Case A, we have z l in the limit [z 0 , z q ] while for Case B, we have z l in the limit [−z 0 , −zq]. Any changes in sign occuring due to change of the order of limits of integration have already been accounted for in writing Eq. H5. The coefficients multiplying the integrals in Eq. H5 can be written as The explicit expressions for these are , C B 2 = z q 2(1 − z γ ) , C A 3 = C B 3 = 2z q q − γ α γ µ , C A;ij 4 = C B;ij 4 = 2 zq γ i γ α γ j γ µ − (1 − zq) γ i γ α γ µ γ j − z q γ α γ i γ µ γ j q − , We can now express the arguments V ⊥ and ∆ of the constituent integrals appearing in Eq. H5 in an expansion in z l , just as in Eq. 154. It is then straightforward to isolate the logarithmically divergent pieces in Eq. H5. As expected, we will only require the expressions for the coefficients c 1⊥ and c 3 that appear in Eq. 154 for the two cases A and B. The rapidity divergent contribution from (V 21) can now be finally written as Recall that what we call the "remainder" (Sβ) is comprised of terms in the divergent part R (Sβ) (I);µα of the amplitude which do not contain UV and ln(1/z 0 ) singularities.
Since the divergent pieces are zero for the processes (S2) and (S3) so are the remainder terms for them. We therefore have (Sβ) µα (l 1⊥ ) = 0, for β = 2, 3 . (I3) For the processes (S1), (S4) and (S6), the remainder has a generic structure that can be written as     , and The coefficients needed for the evaluation of the above terms are given in Eqs. I7 and I8.
We will now present the expressions for the pieces R (Sβ) (II) in the finite part of the amplitude given by Eqs. 102 and 103. These are comprised of integrations over the transverse gluon loop momentum, l 3⊥ , which are finite in two dimensions. We will first write down these contributions for (S1), (S4), (S5) and (S6) which share a similar structure. This will be followed by the corresponding expressions for processes (S2) and (S3) which are considerably more tedious albeit similar to one another.

Vert.(1) finite
The finite pieces for the 6 diagrams in Fig. 14 where R Vert. (1) finite;µα (l 1⊥ ) = In the above equation, for a particular process β, as previously for the self-energy contributions, (V β) is defined as the remainder between terms in the amplitude which are proportional to 1/z l and those proportional to logarithms in z 0 . The second term R (V β) (II) is comprised of finite terms that are not proportional to 1/z l . We will represent the latter in terms of constituent integrals described in Appendix B. For the computation of this second finite piece, it will be useful to express the amplitude for each process (see below) , β = 1, . . . , 6 , (I40) in terms of some generic forms which will be identical for processes with similar topology. We will demonstrate this in detail for a few graphs in the upcoming discussion. The finite pieces for the other graphs can be obtained following similar techniques. These expressions are lengthy and are necessary only in the context of a numerical computation. Mathematica scripts are available that allow for an algorithmic evaluation of these. We begin our discussion with the processes (V 1) and (V 6) which have similar structures with respect to the emission vertex of the outgoing photon, which also extends to their amplitudes. This is exhibited in the similar forms for the rapidity divergent structures in Eqs. 160 and G6.
To compute the first component, (V β) (β = 1, 6) of the finite pieces we need to subtract the divergent part R (V β) div.
from the terms in the amplitude proportional to 1/z l which is represented by R For (V 1), these are given respectively by Eqs. 160 and 151. While performing the integration over z l in the term R (V 1) (I) we will encounter the coefficients c 1 , . . . , c 5 in terms of which the arguments V ⊥ and ∆ of the constituent integrals are expressed as polynomials of z l (see Eq. 154).
In terms of the constituent integrals given in Eq. G10 we can now write the finite pieces contained in R The expressions for theF 's can be obtained in terms of the coefficients b (V β) i (i = 1, . . . , 26) for the processes (V 3) and (V 4). In a similar fashion, one can obtain the finite pieces for (V 2) and (V 5) in Fig. 26. For these, we need to compute only R (II) which also contain the same constituent integrals as in Eq. I57. Similarly to the (V 3) case, we need to include the contribution from l − 3 in the range k − < l − 3 < k − + k − γ in addition to 0 < l − 3 < k − to obtain the net finite contribution.

Vert.(2) finite
The finite pieces of the contribution from the quark-antiquark interchanged counterparts of the six processes in Fig. 26 are contained in M Vert. (2) finite;µα = 2π(eq f g) 2 where R Vert. (2) finite;µα (l 1⊥ ) = We can obtain the net finite contribution from these diagrams by using Eq. 66 along with a change of sign on the various terms constituting the finite piece given by Eq. I39.

Vert.(3) finite
The finite contributions to the amplitudes from the processes (V 13) − (V 15) in Fig. 28 and their q ↔q counterparts are expressed as M Vert. (3) finite;µα = 2π δ(1 − z v tot )(eq f g) 2 dΠ LO ⊥ u(k) In order to obtain these pieces for (V 13) and (V 14) and their q ↔q counterparts, we need to compute the finite pieces in the loop contribution described by the graph in Fig. 35 for the cases l − 2 > 0 (denoted by Case A in the discussion of Sec. V D) and l − 2 < 0 (Case B) and add them up. The computation of the finite contribution from (V 15) is considerably tedious and explicit expressions will only be provided for numerical studies in future. Here we sketch the structure of these pieces for the graph (V 13) in terms of the constituent integrals that appear in virtual graph computations as we have done in the previous sections. The contributions for (V 14) have a similar structure to that of (V 13) and can be obtained following the same methods.
For (V 13), we discussed the computation of the divergent piece in detail in Sec. V D. From the expressions for the gluon loop contribution obtained separately for the two cases in Eqs. 174 and 175, we can isolate the pieces proportional to the constituent integrals I (2,i) v (V ⊥ , ∆) and I (2,0) v (V ⊥ , ∆) which will yield finite results in d = 2 dimensions. A straightforward application of the identity in Eq. D3 in these equations tells us that there are a lot of such contributions proportional to these integrals each with different gamma matrix structures. We can however collect all these pieces and express them as polynomials in z l which will assist in the numerical evaluation of such contributions. With this strategy in mind, we can write the finite contribution to (V 13) as