Inclusive prompt photon-jet correlations as a probe of gluon saturation in electron-nucleus scattering at small $x$

We compute the differential cross-section for inclusive prompt photon$+$quark production in deeply inelastic scattering of electrons off nuclei at small $x$ ($e+A$ DIS) in the framework of the Color Glass Condensate effective field theory. The result is expressed as a convolution of the leading order (in the strong coupling $\alpha_{\mathrm{s}}$) impact factor for the process and universal dipole matrix elements, in the limit of hard photon transverse momentum relative to the nuclear saturation scale $Q_{s,A}(x)$. We perform a numerical study of this process for the kinematics of the Electron-Ion Collider (EIC), exploring in particular the azimuthal angle correlations between the final state photon and quark. We observe a systematic suppression and broadening pattern of the back-to-back peak in the relative azimuthal angle distribution, as the saturation scale is increased by replacing proton targets with gold nuclei. Our results suggest that photon+jet final states in inclusive $e+A$ DIS at high energies are in general a promising channel for exploring gluon saturation that is complementary to inclusive and diffractive dijet production. They also provide a sensitive empirical test of the universality of dipole matrix elements when compared to identical measurements in proton-nucleus collisions. However because photon+jet correlations at small $x$ in EIC kinematics require jet reconstruction at small $k_\perp$, it will be important to study their feasibility relative to photon-hadron correlations.


I. INTRODUCTION
It is conjectured [1,2] that the proliferation of soft (small x) gluons in hadron wavefunctions at high energies via bremsstrahlung is saturated by many-body screening and recombination effects when gluon occupancies are parametrically O(1/α s ), where α s is the QCD coupling. The discovery and characterization of this gluon saturation phenomenon is a major goal of deeply inelastic scattering (DIS) experiments at the Electron-Ion Collider (EIC) and several signatures of gluon saturation have been discussed in this context [3,4].
We will discuss here the process e + A → e + γ + jet + X corresponding to the computation of a photon+quark inclusive final state in DIS off a nucleus with mass number A, where the quark fragments to a jet and an arbitrary number X of other final state hadrons are produced. This process (isolated photon+jet) was previously measured at HERA in DIS off protons [5] and compared to perturbative QCD (pQCD) computations 1 to O(α 2 em α s ) in [8][9][10]. Our focus here is to explore the impact of gluon saturation on photon+jet correlations at the EIC. The distinguishing feature of this phenomenon is an emergent semi-hard saturation scale Q s,A (x) that increases with decreasing x and increasing mass number A. We will study how the saturation scale scale modifies photon+jet correlations. Because all-twist many-body correlations become important in this regime, this requires one to go beyond leading twist perturbative computations at small x and resum powers of Q 2 s,A /Q 2 to all orders at each order in the weak coupling expansion in α s . The physics of "saturated" gluons in QCD is described by the Color Glass Condensate (CGC) effective field theory (EFT) [11][12][13][14][15][16][17][18][19][20] which employs a Born-Oppenheimer separation of the hadron wavefunction whereby the dynamics of small x partons is sourced by static large x parton color charges. The saturation scale Q 2 s (x) appears as a dimensionful scale in this EFT and corresponds to the density of color sources that interact coherently with the small x partons. The density of color sources grows with mass number, thus one obtains Q 2 s ∝ A 1/3 . Since the saturation scale is the only dimensionful scale at x (and A) where Q 2 s (x) Λ 2 QCD , the typical momenta of partons is peaked at this scale. Since then α s (Q 2 s ) 1, the EFT provides a systematic semi-classical weak coupling expansion in the strongly correlated, non-linear small x regime of QCD [11][12][13] where the phase space occupancy of gluons ∼ 1/α s (Q s (x)) 1 is large.
Testing the quantitative reliability of perturbative calculations in the CGC EFT requires one to consider observables which can provide a clear picture of the underlying dynamics. Inclusive prompt photon production and prompt photon-jet angular correlations are two examples of such observables. These processes have been studied thus far in the CGC EFT only in the context of proton-proton/nucleus (p + A) collisions [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37]. In this paper, we will present the first computation of the leading order (LO) cross-section for the inclusive production of a prompt photon in association with a quark (jet) in electron-nucleus (e + A) scattering, in the Regge-Gribov kinematics of fixed Q 2 , squared center-of-mass energy s → ∞ and Bjorken x Bj → 0. Our work employs the formalism for inclusive prompt photon+quark-antiquark (γ +qq) dijet production to next-to-leading order (NLO) developed in [38][39][40]. We will adapt this formalism to study the impact of gluon saturation on inclusive photon+quark correlations. Due to the analytical and numerical complexity of the computations, we will only consider the LO formalism here; the extension to NLO will be left to future work.
At LO in α s , inclusive prompt photon+jet production in e+A scattering at small x proceeds through the fluctuation of the exchanged virtual photon into a long-lived quark-antiquark dipole plus a photon emitted from either parton 2 . The quark and antiquark can multiple scatter coherently off the dense gluons in the nucleus either before or after emitting the photon. As a result of this multiple scattering, the quark and/or antiquark acquires a Wilson line color phase. Indeed, in the eikonal approximation valid at high energies, spatial correlators of these Wilson lines incorporate all the information about the nuclear target accessed in the DIS process at small x. Specifically, in the case of photon+dijet production, the quantum expectation value for the process depends on two-point "dipole" and four-point "quadrupole" correlators [38]. The former is accessed in fully inclusive DIS while extraction of the latter requires more differential measurements such as inclusive dijet production [41,42] and the photon+dijet process discussed in [38].
An interesting question we will address is whether the photon+quark process at LO is sensitive to both dipole and quadrupole correlators. This is not evident from a brute force integration of the phase space of the antiquark in the expression for the cross-section for photon+quark+antiquark production. Towards this end, we will show that the results of [38] can be reexpressed in a form where, for a clearly identified photon+quark hard process, the cross-section is sensitive only to dipole correlators. As a corollary, as the phase space for isolation cut on the photon is reduced, the cross-section will show increasing sensitivity to quadrupole correlators.
In addition to addressing the formal question on the relative sensitivity of dipole and quadrupole correlators, we will explore the discovery potential of this measurement in the kinematics of the EIC. The prospect of such measurements is exciting since they will be performed in e + A collisions at the EIC for the first time. Further, the EIC will operate at unprecedented luminosities which opens up the possibility to measure rare events through comprehensive analyses of experimental data hitherto unavailable even in electron-proton collisions. While identical measurements are feasible in p + A collisions at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC), final state interactions in such experiments complicate the unambiguous extraction of gluon saturation in the hadron wavefunctions. However because the computation of this process also exists in the CGC EFT, comparisons of the results of e+A and p+A collisions offer one the opportunity to systematically isolate initial and final state contributions and thereby extract universal features of gluon saturation. Specifically, we will be able to determine if the dipole and quadrupole operators extracted from different processes in DIS are the same as those extracted in p + A collisions.
Turning now to empirical probes, we note that measurements of azimuthal angle correlations of final state particles are especially sensitive to gluon saturation. In conventional perturbative QCD, energy-momentum conservation dictates that a hard photon (or jet) is produced back-to-back with another jet because only transverse momentum k ⊥ ≤ Λ QCD is assumed to be transferred from a hadron in the computation of the hard cross-section. However, as we noted, with the onset of gluon saturation, the momentum transfer from the hadron is given by k ⊥ ∼ Q s (x). This contributes to a systematic suppression and broadening of the back-to-back or "away side" peak in the ∆φ relative azimuthal distribution of the hard photon-jet cross-section or of dijets.
Much of the discussion on azimuthal angle correlations at small x in the CGC and related frameworks has been restricted to dijet/dihadron production in p(d) + A collisions [41,[43][44][45][46][47][48][49][50][51][52][53][54][55]. Similar studies have been performed for the case of e + A scattering for the production of dijets/dihadrons [42,[56][57][58] and trijets 3 [61]. We will explore here such correlations between the observed photon and quark in the context of e + A collisions at small x. We will show that our results are similar to those for inclusive photon+hadron production in p + A collisions in LHC kinematics [27,28,30,33,62] where a suppression and broadening of the "back-to-back" peak in the azimuthal angle distribution is predicted.
The outline of the paper is as follows. In Section II, we begin with general definitions for the various kinematic quantities involved in this DIS process of interest. We next formulate the computation of the amplitude and elucidate the different components necessary for obtaining the cross-section for inclusive prompt photon+quark production in e + A DIS at small x. The key expression for the cross-section is written in Eq. (18) in terms of contributions from different polarization states of the exchanged virtual photon. The most important ingredients in deriving this quantity are the amplitude and cross-section for γ + qq production in the hadronic subprocess (virtual photon-nucleus scattering), for a given polarization of the virtual photon.
We discuss the computation of the former in section III; the expression for the amplitude is provided in Eqs. (34) and (35). Using this expression for the amplitude, we derive in Section IV, results for the cross-section for γ + qq production separately for the longitudinal and transverse polarizations of the virtual photon. The contribution from the longitudinally polarized photon to the cross-section is given by Eqs. (52) and (54) and the net contribution from the two transverse polarization states is given by Eqs. (59) through (63).
In Section V, divided into three subsections, we obtain results for the cross-section for γ + q production using expressions for the γ + qq cross-section obtained in the previous section and integrating over the phase space of the antiquark. We consider the case of a longitudinally polarized virtual photon to demonstrate this procedure. The corresponding exercise for transverse polarizations is discussed in Appendix F. In Section V A, we discuss in detail the treatment of the singular contribution for the case when the emitted photon becomes collinear to the integrated antiquark. We use a simple cutoff procedure to regulate the collinear divergence and absorb it into the (anti)quark-to-photon fragmentation function; the cross-section proportional to this is given by Eq. (76). In Section V B, we combine the non-divergent terms from the integration over the antiquark phase space into what we call the "direct" photon+quark contribution to the cross-section. The results are given in Eqs. (80)- (85). We then consider in Section V C a limiting form of the "direct" photon+quark production cross-section in momentum space. We show that this contribution is dominated by photons that are hard relative to Q s . This momentum space expression for the direct photon+quark cross-section is given in Eq. (102) with its constituents defined in Eqs. (103)-(107).
One of the main advantages of the simplification in Section V C is that it significantly aids in the numerical analyses for direct photon+quark production. In Section VI we employ these results to study for the first time the azimuthal angle correlations in direct photon-quark in e + A scattering at small x. Since we are able to identify a kinematic region where only dipole Wilson line correlators contribute, we can study the energy evolution of the azimuthal angle correlations by computing the x evolution of this correlator. This is to a good approximation described by the Balitsky-Kovchegov (BK) renormalization group (RG) equation [63,64] that systematically resums leading logarithms containing all powers of α s ln x. In our work, we will employ the running coupling (BK) equation that involves a specific scheme for running coupling effects that appear at next-to-leading log accuracy in the RG evolution [65].
In section VI A we review the setup of the calculation and in section VI B we present the numerical results. We focus our attention on the suppression and broadening of the back-to-back peak as we change the nuclear species from proton to gold (Au). We study the dependence of these observables on different kinematic variables, the transverse momenta and rapidity of of the final state photon and the Q 2 dependence of the virtual photon probe in the scattering process. We summarize our principal results in Section VII and briefly discuss future work in the framework of this computation.
Appendices A-G supplement the material in the main text. The conventions for light cone coordinates are given in Appendix A. Appendix B discusses the Feynman rules used in the CGC EFT calculation. In Appendix C, we present expressions for the squared amplitude for longitudinal and transverse polarizations of the virtual photon that are used to construct the perturbative components of the cross-section in Eqs. (52) and (59). In Appendix D, we list the identities involving traces of gamma matrices and photon polarization vectors that are used to derive the results in the previous section. These are very general and may be useful for other small x computations. Appendix E contains explicit expressions for the various quantities that constitute the expression for the cross-section for γ + qq production from a transversely polarized virtual photon given by Eqs. (59)- (63).
Appendix F details the computation of the cross-section for photon+quark production for the case of transversely polarized virtual photons. Here the discussion is organized in an identical manner to the discussion in Section V and results for the fragmentation and direct photon+quark contributions are presented separately over three subsections. In Appendix F 1, we present in Eq. (F3) the combined contribution from the two transverse polarization states of γ * towards the fragmentation photon+quark cross-section. The corresponding results for the direct photon+quark production are provided in Appendix F 2 in Eqs. (F6) through (F16). In Appendix F 3, we obtain a simplified expression for this cross-section fully in momentum space; this is given by Eq. (F27), with its constituents defined in Eqs. (F28) through (F39). We will use the results of Appendix F 3 in our numerical analyses in Section VI. Finally, Appendix G contains expressions for inclusive quark production in DIS at small x Bj for different polarizations of the virtual photon. These are used in the normalization of the observable that is numerically evaluated and depicted graphically in Section VI.

II. GENERAL CONSIDERATIONS AND OUTLINE OF THE COMPUTATION
In this section, we will discuss features of the CGC EFT framework necessary to compute the cross-section for inclusive prompt photon+quark production in e + A scattering at small x. These considerations are quite general and apply to any inclusive process in e + A collisions.
We begin with a brief description of the relevant kinematics. A relativistic energetic electron with four momentum, k e is scattered off an ultrarelativistic nucleus with momentum P N , per nucleon. The scattered electron has momentum k e and the invariant mass of the recoiling system is given by W . We will work in a frame where the nucleus is right moving and has a large '+' component of nucleon light cone 4 (LC) momentum P + N and both the nucleon and the virtual photon exchanged in the scattering have zero transverse momentum, We assume in this work that P − N = where M N is the mass of the nucleon. The virtuality Q 2 of the virtual photon is completely determined by the four momenta of the incoming and outgoing electron, We will consider the deeply inelastic scattering (DIS) region of e + A collisions, where Q 2 For the kinematics of interest, the largest energy scale in the process is the center-of-mass energy, √ s, where s = (k e + P N ) 2 = 2k − e P + N since we will consider only x Bj = Q 2 sy ≤ 10 −2 , with fixed inelasticity y = P ·q P ·k . We will neglect the mass of the electron throughout the calculation.
As discussed in the introduction, the inclusive prompt photon+quark production in DIS at small x proceeds through the hadronic subprocess in which the exchanged virtual photon (γ * ) fluctuates into a long-lived quark-antiquark (qq) dipole, which multiple scatters off the classical background field of the nucleus, and the quark or antiquark emits a photon (γ) prior or subsequent to this scattering. In order to calculate the cross-section for the production of photon+quark, we need to first compute the cross-section for γ + qq production and then integrate out the phase space of the antiquark.
The 4-momentum assignments for this process are shown in Table I. Following [38], we can write the amplitude for the electron-nucleus scattering process e(k e ) + A(P N ) → e (k e ) + q(k) +q(p) + γ(k γ ) + X , as where X represents anything the nucleus breaks into. Here we have represented respectively by M lep and M had the amplitudes for the leptonic and hadronic subprocesses. The former can be obtained as Π µα (q) denotes the gauge dependent expression for the numerator in the propagator for the virtual photon. We will specify our choice of gauge in section III.
We will now decompose the amplitude in Eq. (3) in a basis of polarization vectors of the virtual photon. One of the advantages of decomposing the amplitude in this new basis is that it allows us to examine the behavior of observables as a function of the polarization of the initial state for the hadronic reaction, i.e., the virtual photon. Although completely equivalent, this is different from the computation in Ref. [38] where the decomposition was in terms of the basis vectors in Minkowski space 5 .
We begin with the well known completeness relation whereby we express the quantity Π µα (q) as a sum over the three polarization states of the virtual photon, We will always use λ (and variants λ ,λ) to indicate the polarization states of the photon (real or virtual).
Taking the modulus squared of the amplitude in Eq. (3), and performing the necessary averaging and sum over electron spins and final state photon polarizations, we obtain where L µν and X αβ are the lepton and hadron tensors, respectively and we have used the completeness relation in Eq. (5) to define in the basis of polarization vectors. One can interpret these as matrix elements of two 3 × 3 matrices L and X , whose diagonal and off-diagonal elements represent the cases of having identical and non-identical polarization states of the virtual photon in the amplitude and the complex conjugate amplitude respectively. The lepton and hadron tensors can be expressed in terms of the amplitudes for their subprocesses as 6 We can therefore rewrite the quantity in Eq. (8) as where is the hadron amplitude for a given polarization state of the virtual photon.
In the CGC EFT, the expectation value of a generic operator O is expressed as where O[ρ A ] is the quantum expectation value for a particular configuration ρ A of the color sources 7 and W Y [ρ A ] is a stochastic gauge invariant weight functional representing the distribution of color charge configurations at a rapidity Y = ln(x 0 /x) in the target 8 . We can therefore write the general expression for the differential cross-section for inclusive photon+quark production in the CGC EFT as 5 In the polarization basis we also have four basis vectors, three of them being the polarization vectors, λ (q), λ = 0, ±1. The fourth one, which is along the four momentum q, gives zero upon acting on both the lepton and hadronic amplitudes by virtue of the Dirac equation and Ward identity, respectively. 6 Here the short form 'pols.' refers to the sum over polarizations of the final state photon. 7 For a general discussion of field theories with strong time dependent sources, see for example [66][67][68] and references therein. 8 x 0 0.01 refers to the initial scale separating sources from fields in the target wavefunction. x < x 0 is the momentum fraction of gluons in the target probed by the projectile and is determined from the kinematics. This is the scale up to which the sources must be evolved to incorporate quantum effects. While it is usually different from the Bjorken variable xps BJ = Q 2 /2P · q, the two can be equated to leading logarithmic accuracy.
where dΩ represents the phase space of the final state particles. For the process of interest, we have chosen to integrate over the phase space of the antiquark (q). For the particles produced from the hadronic subprocess (q,q, γ) we can express their invariant phase space as where in general dη = dv − /v − , with η the rapidity of a particle of four momentum v in the frame we are working in. For the quark, this can be written 9 as η q = 1 2 ln(k − /k + ). Similar relations follow for the antiquark and the emitted photon. In our convention, the rapidity is defined to be positive in the direction of left moving particles, namely, moving along the '−' LC direction 10 .
The CGC averaged (see Eq. (13)) square of the amplitude in Eq. (6) can be written in the basis of polarization vectors as We can express the phase space for the outgoing electron as where φ e is the azimuthal angle of the outgoing electron's 3-momentum, k e . By evaluating the different matrix elements L λλ , using Eqs. (7) and (9), it is easy to show that integrating over φ e projects out only the products of the diagonal elements in the sum over polarizations of γ * appearing in the r.h.s of Eq. (16). From the lepton side, this means that only the elements L 0 0 and L ±1 ±1 enter the expression for the φ e -integrated cross-section.
We can now write the general expression for the differential cross-section for inclusive γ + q production in e + A DIS at small x as where the leptonic contribution is contained in the longitudinally and transversely polarized photon fluxes defined respectively as These are constituted, respectively, of the components L 0 0 and L ±1 ±1 of the lepton tensor defined by Eqs. (7) and (9). The dependence on W 2 enters in the above quantities through the inelasticity variable y we defined previously which can be written (neglecting the mass of the nucleon) as Since we are interested in the kinematics s, W 2 Q 2 , we have y ∼ W 2 /s. In Eq. (18), we have separated the leptonic contribution to the γ + q production cross-section from the hadronic part. The hadronic contribution to the cross-section is contained in the quantity, dσ γ * λ A→qγX , which represents the γ * + A scattering for a given polarization λ = L(0), T (±1) of the virtual photon 11 . In terms of the cross-section for γ + qq production in the virtual photon-nucleus scattering, this can be written as 9 Note that we are working in the massless limit, in which the rapidity and pseudorapidity are identical. 10 We remind the reader that we chose the target nucleus to be right moving as opposed to a convention frequently adopted in the literature where the nucleus is a left mover and the rapidity (of a particle with 4-momentum v) is defined as η = 1 2 ln v + /v − , which is positive along the direction of right moving particles. 11 In what follows we will always sum over the two transverse polarizations T (±1), instead of averaging as typically done.
where dσ γ * λ A→qqγX The CGC averaged hadron tensor can be written in the basis of polarization vectors as where M λ had (see Eq. (12)) is the amplitude for the hadronic subprocess γ * A → qqγX for a given polarization λ of the virtual photon.
We summarize below the steps necessary to obtain the cross-section for γ + q production in e + A DIS at small x: 1. Compute the total amplitude, M λ had for γ + qq production in virtual photon-nucleus scattering. 2. Compute the differential cross-section for γ + qq production as defined in Eq. (22) using the expression for the amplitude.
3. Integrate over the phase space of the antiquark (p ⊥ and ηq) to get the differential cross-section for γ + q production as defined in Eq. (21).
Over the next sections, we will sequentially demonstrate the realization of these steps.
We will work here in the basis of polarization vectors of the virtual photon which was discussed in the previous section. To further simplify the computation, and with a view to extending the current work to higher orders in α s , we choose to work in the light cone gauge A − = 0 using momentum space Feynman rules (see Appendix B).
We choose for the polarization vectors of the virtual photon that we introduced in the previous section, the following convention: It is easy to check that these vectors satisfy the completeness relation in Eq. (5) as λ=0,±1 where we have provided an explicit expression for the numerator Π µα of the photon propagator in A − = 0 gauge. As shown in Fig. 1, at leading order in the CGC power counting, there are four contributions to the amplitude [38], The Feynman graphs along each row in Fig. 1 differ from each other in terms of the photon emission relative to the scattering off the background classical field by the quark-antiquark dipole. Also, the contributions from M 3 and M 4 are obtained from those of M 1 and M 2 , respectively, by interchanging the quark and antiquark lines. The momentum space expression for the amplitude of the hadronic subprocess labeled M 1 can be written as 12 where i,j, are color indices.
Further, q f represents the fractional charge of the quark/antiquark in units of the elementary charge e, the free (massless) quark Feynman propagator is given by and the all-twist resummed effective vertices that enter the momentum space dressed quark and antiquark propagators are represented by Here z ⊥ is shorthand for d 2 z ⊥ andŨ +1 ,Ũ −1 denote lightlike Wilson lines in the fundamental and antifundamental representation of SU (N c ) respectively. Their labels i and j represent color indices in the fundamental representation of SU (N c ) and sign(p − ) is the sign function. The latter takes values of ±1, respectively for the case of quark and antiquark. From the explicit expression provided in Eq. (B4), one can see that the functional dependence on the color charge density ρ A enters the hadron amplitude in Eq. (27) throughŨ . Note that the effective vertices on the quark and antiquark lines are proportional respectively toŨ (x ⊥ ) andŨ † (y ⊥ ), where x ⊥ and y ⊥ are the respective transverse coordinates. The effective vertices given by Eq. (30) represent all possible scatterings of the quark-antiquark dipole with the classical field of the nucleus including the case of "no scattering" [38] (see also Appendix B). Since we are interested in the physical amplitude which should necessarily involve interactions with the nucleus, we have to subtract from the expression for the amplitude in Eq. (28) (and similarly for the amplitudes of the other processes), the possibility of "no scattering" with the background field. This contribution is obtained by settingŨ = 1 in the expressions for T 's that enter the amplitudes for each process in Fig. 1.
Subtracting the "no scattering" contribution from the amplitude in Eq. (28), and using the expressions in Eqs. (29) and (30), we can factor out the color structure and the phase of the quark and the antiquark to write the amplitude for M 1 as where Above we have introduced r ⊥ = x ⊥ − y ⊥ , and the numerator is given by, The integration over l − can be performed trivially to yield an overall momentum conserving delta function: . This is a common feature of all CGC computations performed within the eikonal approximation 13 . We next perform the integration over l + using Cauchy's theorem of residues for complex contour integration. The numerator given by Eq. (33) is independent of l + indicating that the value of the integral over l + will be entirely given by the residue theorem, in the limit of the contour radius going to infinity. The location of the l + poles, in this case, are determined by the signs of the external longitudinal momenta, k − , p − and k − γ . After these steps, we will finally express N λ 1 in terms of a transverse integration over l ⊥ . The amplitudes for the remaining three diagrams can be obtained following the same routine.
We can factor out the overall momentum conserving delta function to rewrite the total amplitude from the four processes in Fig. 1 as where For the process M 1 , we can write where We introduced in this expression the ratios of the outgoing momenta to the component q − of the incoming virtual photon momentum 14 , In deriving Eq. (37) from Eq. (33), we have made use of the Dirac equation:ū(k)/ k = 0, the photon transversality condition, k γ · λ (k γ ) = 0 and the Clifford algebra identity given by Eq. (A2). There are a couple of advantages in writing the amplitude in the form given by Eq. (37). The first one is apparent when one expands the Dirac structure in Eq. (37) as For the choice of polarization vectors in Eqs. (24) and (25), one can easily see that the amplitude naturally lends itself to a decomposition in the basis of polarization vectors. For longitudinal polarization of the virtual photon, only the first term on the r.h.s of Eq. (39) contributes to the amplitude, whereas only the second term survives for transverse polarization of γ * . This feature aids in the computation of the squared amplitude in this basis.
The second advantage of Eq. (37) becomes manifest when we square this amplitude and sum over the polarizations of the final state photon. In this process, one generates from the numerator of the square of Eq. (37) a term |z q k γ⊥ − z γ k ⊥ | 2 that cancels an identical factor in the denominator. This makes the (collinear) singularity for z q k γ⊥ = z γ k ⊥ transparent in the resulting expression.
Through a similar exercise used to derive the amplitude A λ 1 in Eq. (37), we can write the amplitude for the process labeled M 2 in Fig. 1 as where In writing Eq. (40), we have split the amplitude into two contributions labeled 'regular' (reg) and 'instantaneous' (ins). The Dirac structure for the 'regular' term in Eq. (41) is similar to that of A λ 1 (l ⊥ ) in Eq. (37). The nomenclature for the second contribution follows from light-cone perturbation theory, in which the 'instantaneous' term corresponds to the emission of the photon from the qq splitting vertex. This can be seen from the Dirac structure of Eq. (42) where there is no propagator between the two photon vertices.
Because of their similar Dirac structure, the decomposition in Eq. (39) also applies to the 'regular' term in Eq. (41). Using the identity one may verify that the 'instantaneous' contribution is non-zero only for transverse polarizations of the virtual photon. Finally, there is no collinearly singular structure in A 2 , because the photon is emitted by the (off-shell) quark before it scatters off the shock wave.
In this section, we have presented the expressions for the amplitudes of diagrams M 1 and M 2 . For the momentum assigments shown in Fig. 1, the expressions for A λ 3 (x ⊥ , y ⊥ ) and A λ 4 (x ⊥ , y ⊥ ) can be obtained from the corresponding expressions for A λ 1 (x ⊥ , y ⊥ ) and A λ 2 (x ⊥ , y ⊥ ), respectively, by implementing the following replacements in Eqs. (36) and (40): x ⊥ ↔ y ⊥ , p ↔ k (and v(p) ↔ū(k)), an overall '−' sign, and reversing the order of the spinor structures.
In order to obtain the differential cross-section for γ + qq production, we have to multiply the net amplitude from the four processes in Fig. 1, as written in Eq. (34), with the complex conjugate amplitude, and then perform the CGC averaging over the color sources. The seemingly problematic (2π)δ(0 − ) infinite factor that arises from the square (34) can be fixed by using a properly normalized wavepacket description [21,38] instead of a plane wave for the incoming electron. This brings in a factor of 1/2k − e . We can represent this schematically as In writing the cross-section for γ + q production in e + A DIS in Eq. (18), we have absorbed the factor of q − /k − e = y (see Eq. (20)) in the definition of the longitudinally and transversely polarized photon fluxes defined in Eqs. (19). The unpolarized cross-section for inclusive γ + qq production in virtual photon-nucleus collisions, defined by Eq. (23), can therefore be written as where X λλ Y can be written in terms of the reduced amplitudes (see Eq. (35)) for the four diagrams in Fig. 1 as Here and henceforth, counterparts of the coordinates and internal momenta in the Hermitian conjugates of the amplitudes carry prime labels. In the second line of Eq. (46), the first sum is over the helicities of the Dirac spinors and polarizations of the final state photon, while the second sum runs over the four diagrams in Fig. 1 giving a total of 16 contributions at the level of the cross-section. These sums generate traces over the gamma matrices contained in the expressions given by Eqs. (37), (41) (42) and their q ↔q counterparts. Only 6 of these contributions are independent; the remaining 10 can be obtained simply by complex conjugation and q ↔q interchange. These perturbatively computable pieces constitute the process dependent LO coefficient function or in the language of Regge theory, the LO 'impact factor" for γ + qq production. The impact factor for photon+dijet production in e + A DIS at small x was recently computed to NLO [39,40]. Similarly, the color traces over the Wilson lines contained in the squared amplitude and the subsequent CGC averaging give rise to the following combination of dipole and quadrupole correlators: where the dipole and quadrupole Wilson line correlators are defined at a rapidity Y as The CGC averaging inherent in these gauge invariant objects encodes non-trivial correlations between the gauge fields contained in the Wilson lines, thereby representing the coherent many-body gluodynamics characteristic of the saturation regime of QCD. These dipole and quadrupole Wilson line correlators appear in a variety of observables for both e + A and p + A collisions and can be thought of as building blocks of high energy QCD. Although intrinsically non-perturbative, the energy or rapidity evolution of these correlators can be derived using perturbative techniques in the CGC EFT, resulting in the well-known JIMWLK 15 [14][15][16][69][70][71] or BK renormalization group equations. However, one needs to adopt a model to describe the non-perturbative initial distribution of the sources that enters the CGC averaging procedure. A simple and well-known model, which is valid for a large nucleus and moderate x ∼ 0.01, is the McLerran-Venugopalan (MV) model [11][12][13]. In this model, the large x "valence" charges are described as static sources on the light cone, which are Gaussian distributed with where and µ 2 A is the average color charge squared of the sources per color and per unit transverse area of a nucleus with mass number A. Explicit expressions for the dipole and quadrupole many-body correlators are available [19,41,72,73] in the MV model. (See also [74][75][76] for computations of higher point correlators.) We will now provide expressions for the differential cross-section in terms of these "soft" (non-perturbative) and "hard" (perturbative) components explained in the previous paragraphs. For the case in which the virtual photon is longitudinally polarized, we can write the differential cross-section for γ + qq production in γ * L + A scattering very compactly as where d(PS 3 ) = dη q dηqdη γ d 2 k ⊥ d 2 p ⊥ d 2 k γ ⊥ represents the three particle phase space measure and z tot = z q + zq + z γ . We can express the momentum fractions of the final state particles in terms of their rapidities as The coefficient function R qqγ L can be written as We have explicitly shown here the momentum and coordinate dependence of R qqγ L . The terms in each line are organized according to their phases and the labels 1, . . . , 4 denote the respective contributions from the four allowed processes in Fig. 1. The topology of the cut-graphs that represent the different contributions along each line in Eq. (54) is similar. The first (fourth) line, for example, corresponds to the cases in which the final state photon is emitted by the quark (antiquark) in the direct and conjugate amplitudes. Similarly, the second and third lines represent the cross-terms in which the photon is emitted by the quark (antiquark) in the direct amplitude and by the antiquark (quark) in the conjugate amplitude.
The vector functions J i L,p , (p = 1, . . . , 4) are proportional to the light cone wavefunctions for the splitting of γ * L into qqγ; these are distinct because of the different locations of the emission of the photon. We have not shown the momentum dependence in these functions for the sake of brevity. This will be clear from their expressions which are provided below. Finally, the ξ's are coefficients that depend on the momentum fractions; these are defined to be We now provide the expressions for J i L,1 and J i L,2 . The expressions for J i L,3 and J i L,4 are respectively obtained simply by interchange of quark and antiquark: k ↔ p and r ⊥ ↔ −r ⊥ .
In Appendix C, we discuss in detail the steps leading to the expression for R qqγ L in terms of the ξ coefficients and the vector functions J i L . We also provide all the necessary results needed to do this exercise for the transversely polarized virtual photon. This involves explicitly computing the squared amplitudes appearing in the modified hadron tensor in Eq. (46) for the different polarization states of γ * and then organizing the resulting expressions in the symmetric form displayed in Eq. (54). (See also in this regard, Eqs. (61)- (63) in the upcoming discussion.) In a similar fashion, we can obtain the differential cross-section for qqγ production in the case of transversely polarized virtual photons, γ * T . Summing over the two polarization states, we get We separate the perturbative component into three parts as which follows from the decomposition of the amplitude in Eq. (40) into a 'regular' and an 'instantaneous' term. Comparing Eqs. (37) and (41), we immediately see that the 'regular' term has a gamma matrix structure identical to the contribution from the graph M 1 in Fig. 1. We demonstrated using the identity in Eq. (39) that there is indeed a piece in this 'regular' structure that yields non-zero contributions for transversely polarized γ * . We also showed there that the 'instantaneous' term contributes only for transverse polarizations of the virtual photon (see Eqs. (42) and (43)). At the amplitude squared level, we will therefore have direct (reg-reg, ins-ins) and cross (reg-ins or ins-reg) terms. It is therefore evident that there are going to be significantly more contributions to the cross-section in the case of γ * T as compared to γ * L . However, we can carefully organize these terms and express them compactly as shown in the following equations. The first term is given by Comparing the above equation with Eq. (54), we immediately notice a similarity in their structures, albeit with tensor functions, J ij T and different coefficients, ζ and ξ. The explicit expressions for these are provided in Appendix E. Similarly, the second and third terms on the r.h.s of Eq. (60) can be written respectively as and The structures of these terms differ slightly because only the amplitudes of processes M 2 and M 4 in Fig. 1 contain 'instantaneous' contributions. Expressions for all the coefficients, κ and functions, J T ins are provided in Appendix E.

V. INTEGRATING THE ANTIQUARK PHASE SPACE
Now that we have the results for the cross-section for γ + qq production separately for longitudinal and transverse polarizations of the virtual photon, we need to integrate over the phase space of the antiquark to obtain the desired cross-section for γ + q production as given by Eq. (21). We will show here the steps leading to the explicit expressions for this quantity by considering the case of a longitudinally polarized virtual photon. The corresponding computation for the case of transverse polarization is considerably more tedious and expressions for the γ * T A → qγX cross-section are provided in Appendix F.

A. Collinear singularity and the fragmentation contribution
To begin with, the integration over the antiquark rapidity (dηq = dzq/zq) in Eq. (21) can be done trivially using the delta function that enters the expression for the cross-section, in Eq. (52). This sets zq = 1 − z q − z γ . For the transverse momentum integration, we notice that all such terms in the γ * L A → qqγX cross-section (see Eqs. (52) and (54)) that do not have an explicit p ⊥ -dependence will yield a two dimensional Dirac delta upon integration over the transverse phase space of the antiquark. In terms of the vector functions constituting the expression in Eq. (54), these would be the ones not proportional to J i L,3 (or its complex conjugate) which we write below for reference, For the remaining terms in Eq. (54), we have those proportional to J i L,3 (or its complex conjugate) and one term proportional to J i L,3 J i * L,3 . We get finite results for the former after performing the integration over p ⊥ . In this subsection, we will mainly be concerned with the latter which contains a singular contribution after the p ⊥ -integration. Physically, the term proportional to J i L,3 J i * L,3 represents the case in which the photon is emitted from the antiquark after the scattering off the nucleus in the direct amplitude, and by the scattered antiquark in the conjugate amplitude. Using the expression for J i L,3 given in Eq. (64), we see that the square of this quantity gives where we have used the identity When we are integrating over p ⊥ , we will encounter a region where the photon is collinear to the outgoing antiquark. This is apparent from the integral Above we have a p ⊥ -dependent phase that enters the cross-section as given by Eq. (52) and the denominator is taken from the squared expression in Eq. (65). There are many ways to regulate the collinear singularity. The simplest way is to introduce an infrared (IR) regulator Λ in order to write the integral over the transverse momentum of the antiquark as 16 Note that 1/|y ⊥ − y ⊥ | acts as a UV regulator for the p ⊥ integration; these transverse coordinates are also present in the fundamental Wilson lines constituting the LO color structure in Eq. (47). The logarithm can be broken into two parts as by introducing a factorization scale µ F . The contributions from momenta below µ F , contained in the second logarithm on the r.h.s of Eq. (69), is absorbed into the (anti)quark-to-photon fragmentation function Dq →γ (λ γ/q , µ F ). Here λ γ/q (0 ≤ λ γ/q ≤ 1) represents the fraction of momentum of the antiquark carried by a photon. In terms of the partonic variables used in our computation, this is given by We will call this contribution proportional to Dq →γ (λ γ/q , µ F ) the "fragmentation" (F) contribution. The fragmentation function Dq →γ (λ γ/q , µ F ) accompanies the cross-section for the corresponding γ * A → qqγX hadronic subprocess in which the antiquark hadronizes and then emits a "decay" or fragmentation photon. We define the (anti)quark-tophoton fragmentation function 17 at µ F as [7] Dq →γ (λ γ/q , µ F ) = α em q 2 where µ 0 ∼ Λ ∼ Λ QCD and the (anti)quark-to-photon fragmentation function at the initial scale, Dq →γ (λ γ/q , µ 0 ) is obtained from experimental data [77]. The familiar leading order photon splitting function is One can think of this as an RG procedure in which the (anti)quark-to-photon fragmentation function at some initial scale, µ 0 is matched on to that at an evolved scale µ F . The quantum corrections are logarithmic in nature and the hard coefficients are given by the splitting functions. 3 proportional to the first logarithm on the r.h.s of Eq. (69) will be absorbed into the remaining terms that are finite after the p ⊥ integration and together they constitute the "direct" (D) contribution 18 . The inclusive prompt photon+quark production cross-section given by Eq. (21) can therefore be expressed as We will discuss the "direct" contribution in the next section. Here we obtain results for the "fragmentation" contribution to the γ + q production cross-section. For this we will use the expression for the γ + qq production cross-section for longitudinally polarized virtual photon given by Eqs. (52) and (54) and the results in Eqs. (68) and (69). It is easy to check that the prefactor ξqq in the term proportional to J i L,3 J i * L,3 , defined in Eq. (55) is proportional to the photon splitting function defined in Eq. (73), For the longitudinally polarized virtual photon, the fragmentation contribution to the inclusive γ + q production can therefore be written as This can be further expressed in terms of the lightcone wavefunction for the longitudinally polarized γ * as 17 There is a theoretical uncertainty related to the choice of the factorization scale µ F . For example, one can choose µ F to be the invariant mass of the photon-antiquark pair or the transverse momentum of the photon. Although this choice is arbitrary, a knowledge of the fragmentation function at a certain scale can be used to evaluate its value at another scale, thanks to the physical cross-section being independent of µ F . This leads to the evolution equation which also has higher order corrections in the coupling, that are not written above. 18 They therefore have a logarithmic sensitivity to µ F that goes away when this contribution is combined with the fragmentation contribution to obtain the physical cross-sectiion.
where we have used the relations in Eq. (53) to express z q in terms of momentum and rapidity of the quark, and the lightcone wavefunction is defined to be [41] Ψ L α β (q − , z q , r ⊥ ) = 2π where z is the momentum fraction carried by the quark from a longitudinally polarized virtual photon with momentum q that splits into a qq dipole with size r ⊥ . Here we have ∆ = Q 2 z(1 − z) for massless quarks and α, β are the helicities of the quark and antiquark. We obtain similar expressions for the case of transverse polarization of the virtual photon; the computation is described in Appendix F 1. We immediately notice a similarity between Eq. (77) and the LO cross-section for inclusive dijet production in γ * L + A scattering [41]. Interestingly, both these cross-sections contain the quadrupole Wilson line correlator, which is sensitive to the Weizsäcker-Williams (WW) unintegrated gluon distribution in the so-called "correlation" limit 19 of back-to-back hard dijets with small momentum imbalance [41,78]. This makes the study of the fragmentation photon contribution for inclusive prompt photon plus jet production in DIS at small x appealing in its own right.
The numerical analysis, however, is more challenging when one encounters the quadrupole. In the context of DIS, numerical results that incorporate the quadrupole (in addition to the dipole) have been recently obtained for inclusive dijet production in Ref. [42]. This means that it is indeed feasible to provide numerical estimates of the fragmentation contribution for inclusive photon plus jet production, albeit with scale and scheme uncertainties associated with the fragmentation function Dq →γ (λ γ/q , µ F ). We will explore this calculation in the future; in this present work, we will limit ourselves to computing the direct contribution in Eq. (74).
To quantify the full impact of the direct contribution, the numerical results that we will present in this paper will be for kinematics in which the fragmentation contribution is suppressed, which as evident from Eq. (77) is the case for large k γ⊥ /λ γ/q .

B. Obtaining the direct photon+quark contribution
In this section, we will derive an expression for the direct contribution to the differential cross-section in Eq. (74). In order to obtain this for the case of the longitudinally polarized virtual photon, we will compute the antiquark phase space integrated (finite) result for the terms in Eq. (52) modulo the piece of d that has been absorbed into the fragmentation function. As discussed before, most of the terms are independent of p ⊥ and hence will lead to y ⊥ = y ⊥ after integration over p ⊥ and y ⊥ . For such terms, the color structure of the Wilson line correlators simplifies to which only contains dipole correlators. We have chosen to denote the coincidence limit y ⊥ = y ⊥ of the LO color structure in Eq. (47) by Θ to avoid confusion with the color structure for the remaining terms, which contain at least one J i L,3 (or its complex conjugate) and therefore possess the most general color structure with both dipole and quadrupole correlators.
The direct photon+quark contribution to the cross-section can therefore be split into the sum of two terms as where A L has the color structure given by Eq. (79), representative of terms in the amplitude independent of p ⊥ , and B L , which has the general color structure given by Eq. (47). The first piece can be expressed as The expressions for the various coefficients ξ qq , ξqq and ξ qq are given respectively in Eqs. (55) and (56). We should equate zq = 1 − z q − z γ in these expressions because we are integrating over the antiquark phase space. The same argument also holds for the functions J i L,1 , J i L,2 and J i L,4 , given respectively by Eq. (57), Eq. (58) and the q ↔q interchanged counterpart of Eq. (58).
It should be noted that if we assume translational invariance in the transverse direction, then for y ⊥ = y ⊥ , the dipoles appearing in the color structure given by Eq. (79) are dependent only on the relative separations, r ⊥ = x ⊥ −y ⊥ , r ⊥ = x ⊥ − y ⊥ , and r ⊥ − r ⊥ = x ⊥ − x ⊥ . We will use this in the next section where we obtain limiting expressions for the direct photon+quark contribution and express them entirely in momentum space.
To obtain a compact expression for the second contribution B L in Eq. (80), we will first isolate the piece in the r.h.s of Eq. (68) that contains ln(1/|y ⊥ − y ⊥ |µ F ). As discussed previously, this comes from the integration in transverse momenta over p ⊥ of the quantity J i L,3 J i * L,3 , for momenta above the fragmentation scale µ F . We will also use the following result 20 for terms in the amplitude squared that contain a single J i L,3 or its complex conjugate, where r yy ⊥ = y ⊥ − y ⊥ and we have introduced the shorthand notation The corresponding expression for J i * L,3 is obtained by replacingĴ L,3 (r ⊥ ) on the r.h.s. of Eq. (82) byĴ * L,3 (r ⊥ ), wherê As before, r ⊥ = x ⊥ − y ⊥ and r ⊥ = x ⊥ − y ⊥ . We obtain finally, TheĴ L,3 's are defined in Eqs. (83) and (84). It is assumed that the replacement zq = 1 − z q − z γ has been made in the expressions for the different ξ's and J 's appearing in the above equation.
The direct photon+quark production cross-section for the case of transverse polarization of the virtual photon can be similarly broken up into two pieces depending on the color structure. Explicit expressions for these are provided in Appendix F 2.

C. Simplifying the direct photon+quark contribution
We will now consider the direct photon+quark production cross-section obtained in the previous section in a limit where the scale 1/(y ⊥ − y ⊥ ) 2 is much larger than the fragmentation scale µ 2 F . In this limit, we are in a region entirely dominated by direct photons. A strong motivation to study this limit is that the color structure of the B L term in the direct photon+quark cross-section (see Eqs. (80) and (85)) simplifies greatly to include only dipole correlators when y ⊥ → y ⊥ . This will allow us to express the direct photon+quark cross-section entirely in momentum space, which makes the result much more amenable to numerical analysis compared to the general result, which contains 20 To avoid confusion, we note that J i L,3 (r ⊥ ) contains an implicit dependence on p ⊥ . the quadrupole correlator. We will again restrict our discussion here to longitudinally polarized virtual photons. The same techniques can be applied to obtain results for transverse polarizations; these are provided in Appendix F 3. As is evident from the constituent expressions for the direct photon+quark production cross-section, we only need to examine the limiting form for the second contribution B L , since A L (in Eqs. (80) and (81)) already contains only dipole correlators. Taking the naive limit y ⊥ → y ⊥ in Eq. (85) is problematic because the logarithm in r yy ⊥ will diverge. We will therefore perform an approximation in which we Taylor expand the color structure around r yy ⊥ = 0 and perform the integration over r yy ⊥ using known identities.
From the expression for B L given by Eq. (85), we see that the integrals of interest can be expressed schematically in the following form where, for our computation, we have k γ⊥ and f (u ⊥ ) represents the color structure denoted by Ξ(x ⊥ , y ⊥ ; y ⊥ , x ⊥ |Y ) in Eq. 47. We can now do a Taylor expansion for the function f around u ⊥ = 0 ⊥ and rewrite the integral in Eq. (86) as We will now use the well-known identity (see for example Appendix A of [19]) where J 0 is a Bessel function of the first kind. In addition, we make the argument that the physical scale that sets the scale for the gradients of the Wilson line correlators contained in the function f , is the saturation momentum scale Q s , i.e. ∂ k f (0 ⊥ ) ∼ Q s . We can then write the integrals in Eqs. (86) and (87) as where the relative magnitude between consecutive terms is O(Q s /v ⊥ ). In the limit, we can therefore consider only the leading term in the expansion and study qualitative features of the cross-section which is proportional to f (0 ⊥ ) or equivalently only contains dipoles as in Eq. (79). We will now use the above arguments to obtain the leading term in the limit given by Eq. (92) for the term B L in Eq. (85). We consider the dipoles to be translationally invariant so that for y ⊥ − y ⊥ = 0, and we can write Finally, we make the following simple redefinitions in order to obtain a compact expression for the direct photon+quark cross-section that may be expressed entirely in momentum space: The functions J i L,1 and J i L,2 are given by Eqs. (57) and (58). J i L,4 is obtained from Eq. (58) by q ↔q interchange whileĴ L,3 is given in Eq. (83). The expressions for the complex conjugates can be obtained straightforwardly from the above equations. In Eqs. (95)-(98) appearing above, K L 's represent the momentum space versions of their coordinate space counterparts.
In terms of these newly defined functions K i L,p (p = 1, . . . , 4), we can combine the terms A L and B L in Eqs. (81) and (85) to obtain the following expression for the leading contribution to the direct photon+quark cross-section-in the limit described by Eq. (92): where S ⊥ represents the transverse area of the target nucleus and We can now define the CGC averaged dipole correlators as and use Eqs. (95)-(98) to express Eq. (99) in terms of momentum space quantities as where In Eq. (102), C Y (l ⊥ ) contains the non-perturbative input that has to be described by a model such as the MV or Golec-Biernat Wusthoff (GBW) [79] models. The MV model provides initial conditions for the nucleus in the CGC approach whereas the GBW model describes high energy DIS data using a simple parametrization for the dipole amplitude.
For the longitudinally polarized virtual photon, the functions K L constituting the momentum space expression for the cross-section in Eq. (102) are obtained as The coefficients ξ's are given in Eqs. (55)-(56) and Similar expressions for transversely polarized virtual photons are provided in Appendix F 3. Substituting these results into Eq. (18), we obtain a simple momentum space expression for the differential cross-section for prompt "direct" photon+quark production in e + A DIS at small x. Using suitable models for the non-perturbative content, we will now evaluate these expressions numerically and study observables that provide novel ways to study gluon saturation.

VI. AZIMUTHAL ANGLE CORRELATIONS IN DIRECT PHOTON+QUARK PRODUCTION
Azimuthal angle correlations in two-particle production have been proposed previously to be a sensitive probe of gluon saturation inside nuclear matter [43]. Gluon saturation leads to a suppression of the back-to-back peak in the correlation of dihadrons/dijets in e + A [56] and p + A [45,52,53,80] collisions, and in the correlation of photons and hadrons in p + A collisions [36,62]. To go beyond these studies, we compute azimuthal angle correlations of photons and quarks in e + A DIS for the first time. We focus our attention on the production of direct photons, and leave the analysis of fragmentation photons for future work. Towards this end, we will study the behavior of the leading contribution to the differential cross-section (normalized by the transverse area of the target) as a function of the azimuthal angle between the direct photon and the quark ∆φ qγ , and integrated over transverse momenta and rapidities: The differential cross-section for the leading contribution to direct photon+quark production can be found by computing Eq. (18) together with Eq. (102) and Eq. (F27). We further normalize our results by the single inclusive quark production cross-section (see Appendix G) to define the correlation coefficient (associated yield) of direct photons for a quark/jet trigger 21 : where we have summed over the contributions of the three light-quark flavors in defining Eqs. (108) and (109).
A. Set-up: proton and nuclear dipole correlators from running coupling Balitsky-Kovchegov evolution Azimuthal angle correlations from the direct photon+quark alone depend on C Y (l ⊥ ), the two-point dipole correlator in momentum space. We obtain C Y (l ⊥ ) from the Fourier transform of the coordinate space solution to the running coupling Balitsky-Kovchegov (rcBK) equation [63,64,[81][82][83]: where r 2⊥ = r ⊥ − r 1⊥ , and the running coupling kernel K rc is given by [84], The coordinate space coupling in the above is given by with Λ QCD = 0.241 GeV, N f = 3, and C a dimensionless parameter controlling the running of the coupling. The initial conditions for the evolution are given by the modified MV model with a regulator e c (MVe) used in prior phenomenological studies of data from HERA and from proton-nucleus collisions at RHIC aand the LHC [85]: The parameters are C = 7.2, e c = 51.4 and Q 2 sp,0 = 0.060 GeV 2 . Note that the latter two correspond to a saturation scale Q 2 s = 0.238 GeV 2 , defined by S 2 Y =0 (r 2 = 2/Q 2 s ) = exp(−1/2). They have been choosen to fit inclusive DIS HERA e+p data, and provide a good phenomenological description of the proton dipole correlator [85]. (See also [65].) The nuclear dipole correlator is computed via the optical Glauber model, generalizing Eq. (114) by the replacement, where A is the total number of nucleons in the nuclear target, S p,⊥ = 16.4 mb is the effective area of the proton 22 , and is the nuclear thickness function obtained by integrating the Woods-Saxon distribution along the z−direction, with n chosen such that b ⊥ T A (b ⊥ ) = 1. For the gold nucleus (A = 197), we take R A = 6.37 fm, and a = 0.54 fm.
Determining the impact parameter in electron-nucleus collisions is difficult due to the lack of centrality classes. We will here, in this exploratory study, evaluate the saturation scale at the mean impact parameter ( b ⊥ = 0.59 × R A ) as a proxy for minimum bias. This results in the relation Q 2 sA,0 ( b ⊥ ) = 2.7 × Q 2 sp,0 between the nuclear and proton saturation scales at the initial scale x 0 corresponding to the initial rapidity Y = 0. The dipole correlator in momentum space can be computed with the given input parameters for both the proton and the gold nucleus. The result is plotted in Fig. 2.
Given these initial conditions, the rcBK equation is evolved up to Y = log(x 0 /x g ), where x g is determined by the kinematics of the final state particles. For q + γ production we have, where E p is the energy of the proton. This relation follows from energy-momentum conservation 23 . We will perform our study in the frame where the virtual photon and the proton have zero transverse momentum and the momentum of the proton is that of the laboratory frame.
B. Numerical results for azimuthal angle correlations in e + p and e + Au collisions We present in Fig. 3 our numerical results for the direct photon-jet azimuthal angle correlations in e+p and minimum bias e + Au collisions at the center-of-mass energy √ s = 90 GeV which is close to the anticipated EIC center-of-mass energy per nucleon. The kinematical range in transverse momenta and rapidities of the photon and quark have been choosen so that x g remains smaller than 10 −2 . Further, we select an inelasticity y close to 1 so that we maximize the center of mass energy squared W 2 of the photon-nucleon system for a fixed s. The left panel of Fig. 3, shows the yield of direct photons associated to DIS events with at least one quark/jet as a function of the relative azimuthal angle between the direct photon and the jet. We observe a back-to-back peak in the azimuthal angle correlation (∆φ qγ ∼ π), which is suppressed and broadened for collisions with the denser nuclear target compared to those in proton. The origin of the supression in the back-to-back peak can be traced to the broadening of the dipole correlator in momentum space as the saturation scale is increased, as shown in Fig. 2. Physically, the broadening of C Y (l ⊥ ) corresponds to an increase in the momentum transfer imparted to qγ system as it multiple scatters from the nuclear target, resulting in the decorrelation of the back-to-back peak. An additional imbalance of the transverse momenta of the γq is produced by the transverse momentum carried by the integrated antiquarkq; however, the integration over its phase-space is dominated by small transverse momenta.
We also observe an enhacement in the correlation as the photon is emitted collinearly from the quark (∆φ qγ ∼ 0, 2π). The collinear enhancement can be easily identified in Eqs. (104) and (F34). It occurs in the contribution where the photon is emitted from the quark after it scatters from the nuclear target. Our numerical result shows that the collinear region is less sensitive to nuclear effects, as little difference is observed when the nuclear species is changed from proton to gold.
To more efficiently characterize the supression in the back-to-back peak, we study the ratio shown in the right panel of Fig. 3. Near the back-to-back peak we observe a significant supression for the yield of direct photons when comparing proton DIS to nuclear DIS in the same transverse momenta and rapidity bins. We can also study the dependence of R eA on the virtuality Q 2 of the virtual photon. We find that the ratio R eA converges towards unity as the Q 2 of the DIS probe increases, signaling the weakening of the suppression of back-toback correlations. This behavior is characteristic of gluon saturation: as the virtuality Q 2 is increased, the system becomes more dilute and the multiple scattering and shadowing of gluons becomes less dominant. In the dilute limit of Q 2 Q 2 s , one can approximate the Wilson line correlators in Eq. (94) by where S In this approximation, the saturation scale appears as a multiplicative factor in the differential cross-sections in Eqs. (99), (F20), (G4) and (G5). It therefore cancels in the correlator C(∆φ qγ ), resulting in the ratio R eA becoming identical to unity. Thus deviations from unity in R eA provide a good measure of the sensitivity of the azimuthal correlations to the many-body multiple scattering and color screening ("shadowing") effects due the dense gluon system comprising nuclear matter at high energies. The study of the Q 2 dependence of correlation measurements provides an additional handle in e + A collions relative to p + A collisions for systematic studies of gluon saturation phenomena. More detailed information can be obtained by studying the azimuthal correlations at different transverse momenta and rapidities of the photon. The panel on the left in Fig. 5 shows the R eA parameter defined in Eq. (118) as a function of the azimuthal angular separation ∆φ qγ for different ranges of the photon transverse momentum in the same rapidity ranges for quark and photon. We find that the suppression is maximal when the transverse momenta of the quark and photon are closer in magnitude or equivalently, when the back-to-back imbalance in transverse momenta is close to zero. As the difference in the magnitudes of the transverse momenta of the observed final state particles is increased, the suppression factor between that for the proton and for the gold nucleus is reduced. Thus precise reconstruction and identification of the jet and photon transverse momenta respectively are needed to carefully measure the back-to-back suppression. In the case of dihadron production in p + A collisions, it has been argued [52] that the rapidity dependence on the suppression provides an additional handle to probe gluon saturation. It is anticipated that because of small x evolution, particle production at forward rapidities will result in more visible effects of saturation. However, we find very little dependence of the suppression on the rapidity of the photon (right panel in Fig. 5). This might be because of the limited phase-space allowed by the kinematics at EIC energies, which results in limited small x evolution. It may also be due in part to the integration over the antiquark phase space which results in a more complicated dependence of the impact factor on the rapidity of the jet and the photon.

VII. SUMMARY AND OUTLOOK
We derived in this paper analytic expressions for inclusive prompt photon+quark production in e+A DIS at small x and at LO in the CGC EFT power counting. We presented a complete calculation accounting for both fragmentation and direct photon contributions. We find that the former is a convolution of the differential cross-section for inclusive dijet production and the antiquark-to-photon fragmentation function; it is sensitive to both dipole and quadrupole Wilson line correlators. In contrast, at large transverse momenta (which effectively impose an isolation cut on the photon), the direct photon contribution only depends on the dipole Wilson line correlator.
We expressed our results for direct photon+jet production fully in momentum space, and computed the dipole correlator using the solution of the rcBK small x evolution equation. This allowed us to study for the first time the inclusive direct photon-jet azimuthal angle correlations in the kinematic range of the Electron-Ion Collider. We observed a significant suppression and broadening of the back-to-back peak as the saturation scale is increased by comparing e + p collisions to e + Au min-bias collisions. Our results suggest that the back-to-back suppression will be further amplified if it were possible to identify more central e + A collisions [86] corresponding to larger values of the saturation scale. While our results suggest that photon-jet correlations are a promising channel to probe gluon saturation, this will require jet reconstruction for k ⊥ values in the range of 2 − 3 GeV at EIC energies.
Since such jet measurements are extremely challenging, photon-hadron correlations (whereby the final state quark fragments into a hadron), may provide an alternative channel that is sensitive to the process we have computed at the parton level. A useful study, outside the scope of the present work, would be to estimate the relative uncertainties in low k t jet reconstruction and those arising from our anticipated knowledge of fragmentation functions in the EIC era.
In future work, we plan to extend our numerical results to NLO in the CGC power counting. This will require the numerical evalution of the impact factor [39,40] and the numerical implementation of the small x evolution of Wilson line correlators employing renormalization group evolution via the NLO BK [82,[87][88][89] or NLO JIMWLK [90,91] equations.
An important subset of NLO contributions are those of Sudakov type that have been studied in the context of two particle correlations at small x [92,93]. They provide an additional suppression of the back-to-back peak in the production of di-hadrons [56] which may complicate the extraction of a clean signal for gluon saturation. We expect that due to the integration over the phase space of the antiquark in our process, Sudakov effects will be realized differently for photon+quark correlations. If so, a comparative study of the dijet and photon+jet channels can help isolate the relative contributions of Sudakov and gluon saturation effects.
An interesting extension of the computation discussed here would be to integrate over the phase space of the quark and extract the inclusive prompt photon cross-section in e + A DIS at small x. This has been computed at small x for p + A collisions at LO [21] and NLO [31,32,35] in the CGC power counting. Because one has clean initial and final states, essentially free from hadronization uncertainties, the measurement of prompt photons at the EIC provides another important channel to probe the gluon saturation regime. We will explore this possibility in the future.
We work in light-cone coordinates suitable for computations in the eikonal limit. These are defined as with transverse coordinates remaining the same as in Minkowski space. We write four vectors as a µ = (a + , a − , a ⊥ ), where a ⊥ is the two-dimensional vector of transverse components of a. The magnitude of the two dimensional vector will always be denoted by a ⊥ . The metric has the following non-zero entries: g +− = g −+ = 1 and g ij = −δ ij (i, j = 1, 2). The same definition holds for γ + and γ − , with the Clifford algebra defined by the anticommutation relation The scalar product is given by We will use the following shorthand notations for the momentum integrals,

Appendix B: CGC Feynman rules
We work in the light cone gauge A − = 0. In addition to the regular Feynman rules that follow from the QED and QCD Lagrangians (see for example, the book by Schwartz [94]), we also need an expression for the dressed quark propagator in the background of the classical field of small x gluons inside the dense nucleus. In our computation, the Feynman diagram for the dressed quark propagator is shown in Fig. 6. The momentum space expression for the dressed propagator can be written as [13,95,96], where we have considered the two fermion lines on either side of the effective interaction vertex, shown by a crossed circle in Fig. 6, to be internal (off-shell). i, j, k, m represent color indices in the fundamental representation of SU (N c ).
For an outgoing on-shell quark line with momentum p in Fig. 6, we will replace S 0 (p) by the Dirac spinorū(p) in Eq. B1. In Eq. (B1), the free massless quark propagator with momentum p is given by The effective vertex appearing in the dressed quark propagator has the following expression in A − = 0 gauge [95,96] HereŨ is a Wilson line in the fundamental representation of SU (N c ) which can be written as where P − means path ordering along the '−' LC direction, t a (a = 1, . . . , N c ) represents the generators of the fundamental representation of SU (N c ) and A + cl is a solution (see Eq. (B6)) of the classical Yang-Mills equations in the A − = 0 gauge. The latter can be written as where the covariant derivative D µ = ∂ µ − igA µ , F µν is the field strength tensor, 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 right moving at nearly the speed of light. The solutions of the YM equations in the light cone (LC) A − = 0 gauge are given by where Λ is an infrared cutoff necessary to invert the Laplace equation −∇ 2 ⊥ A + cl = gρ A . In Eq. (B1), we include the possibility of "no scattering" (given byŨ = 1) within the definition of the effective vertex . So the dressed propagator also contains a free part given by (2π) 4 δ (4) (p − q)S 0 ij (p) and an interacting part which contains all possible scattering with the nuclear shock wave. This is pictorially depicted by Fig. 7. From this equivalent representation and from the expressions for the fundamental Wilson line in Eq. (B4) and the small x classical gluon field given by Eq. (B6), we observe that the dressed quark (and antiquark) propagators in the CGC EFT efficiently resum all higher twist contributions ρ A ∇ 2 ⊥ → Q 2 S Q 2 from the multiple scattering of the qq dipole off the color field of the nucleus.
Appendix C: Derivation of R qqγ L and R qqγ T In obtaining the differential cross-section for γ + q production in e + A DIS at small x, two essential elements are the "hard" coefficient functions R qqγ L and R qqγ T which enter the γ + qq cross-section in Eqs. (52) and (59) respectively. These functions are extracted from the modified hadron tensor given by Eq. (46). Using this definition for the hadron tensor and comparing it with Eqs. (45) and (52), (59), we can write where A λ m (m = 1, . . . , 4) are the reduced amplitudes for the four subprocesses. Their expressions for the processes labeled M 1 and M 2 in Fig. 1 are given by Eqs. (36) and (40). For their q ↔q interchanged counterparts, we simply replace k ↔ p and x ⊥ ↔ y ⊥ in these expressions and put an overall minus sign. Also, we use the prime symbol to label transverse coordinates and momenta in the complex conjugate of these amplitudes.
Firstly, we note that we do not need to evaluate all the combinations of the squared amplitudes because most of them are related by charge and complex conjugation. We can write the reduced amplitudes in Eqs. (36) and (40) schematically as where the denominators D 1 and D 2 (and their quark-antiquark interchanged counterparts) are not important for the current discussion. The terms of interest are the A λ (l ⊥ )'s which contain the Dirac gamma matrix structures specific to the hadronic subprocesses. When we sum over the spins of the outgoing quark and antiquark at the level of amplitude squared, this results in traces over these gamma matrices. In addition, we have to sum over the two transverse polarizations of the outgoing photon which is denoted by the indexλ in our amplitudes-see Eqs. (37), (41) and (42). For the case of the longitudinally polarized virtual photon (λ = L = 0), the computation of the trace can be simplified using the identity in Eq. (39) and also by noting that there is no "instantaneous" contribution from the amplitude A L 2 . We will provide below the final expressions obtained for the six independent terms in the sum appearing in Eq. (C1).
a. Longitudinally polarized virtual photon spins , pols.(λ) spins , pols.(λ) spins , pols.(λ) spins , pols.(λ) spins , pols.(λ) where the vector functions I i 1,2 are shorthand for The expressions for I i 3 and I i 4 are obtained respectively from I i 1 and I i 2 by the replacements z q ↔ zq and k ⊥ → p ⊥ . It is easy to check that the expressions for the remaining combinations of squared amplitudes can be evaluated from the quantities in Eqs. (C4)-(C9).
To obtain the desired form for R qqγ L as it appears in Eq.(54), we have combined the vector functions in Eqs. (C10) (and their q ↔q counterparts) with the denominators D 1 and D 2 in the transverse integration over l ⊥ to define the new vector functions which we call J i L,p (p = 1, . . . , 4). The coefficients ξ's are formed by the functions of z q and zq that survive after cancellation with the prefactors appearing in Eq. (C1).

b. Transversely polarized virtual photon
In the case of the transversely polarized virtual photon, we have additional contributions due to the non-vanishing of the instantaneous term in Eq. (40). In addition to the sums over the four allowed processes, outgoing quark-antiquark spins and polarizations of the produced photon, we have to add the contributions from the two transverse polarization states of γ * to obtain the desired expression for R qqγ T . We will represent this by the shorthand notation, .

(C11)
In the previous section, we provided explicit results for the independent terms in the squared amplitude for both longitudinal and transverse polarizations of the virtual photon. While the computation for the longitudinal case is very straightforward, the same cannot be said for transversely polarized virtual photon. We will present here a detailed list of the traces over various non-trivial structures involving transverse gamma matrices that appear in this computation. We have explicitly indicated the sum over the two polarization statesλ = ±1 of the outgoing photon.
Four gamma matrices between polarization vector structures: Four gamma matrices (two and two) between polarization vector structures: discussed in the main text for longitudinally polarized virtual photon albeit with more contributions in each case.

Fragmentation photon+quark production cross-section
We begin with the expression for the cross-section for γ + qq production in γ * T + A scattering given by Eqs. (59) and (60) and then integrate over p ⊥ and ηq (or zq). The integration over the latter is trivial using the Dirac delta function and sets zq = 1 − z q − z γ . For the integration over p ⊥ , we have mostly terms in R qqγ T which are independent of p ⊥ and therefore yields the delta: δ (2) (y ⊥ − y ⊥ ). Integrating over y ⊥ , we get the reduced color structure containing only dipoles defined by Eq. (79) for such terms.
There are two kinds of terms that depend on p ⊥ : those proportional to the tensor function J ij * T,3 (or its complex conjugate) which yield finite results upon the p ⊥ integration, and one term 24 proportional to ζqq δ im δ jn J ij T,3 J mn * T,3 which will exhibit collinear singularity. Both these terms have the general LO color structure in Eq. (47). We use an IR regulator Λ and introduce a scale, µ F to isolate the contribution from such fragmentation photons. The contributions above the scale µ F are absorbed into the remaining finite terms and together they constitute the direct photon+quark contribution.
In terms of the momentum fraction, λ γ/q of the photon with respect to the parent antiquark, defined in Eq. (70) we obtain where the photon splitting function is defined by Eq. (73). Also, using the identity where K 1 is a modified Bessel of the second kind, we can express the fragmentation contribution to γ + q production in γ * T + A scattering as Dq →γ (λ γ/q , µ F ) λ γ/q z 2 q + (1 − z q ) 2 r ⊥ .r ⊥ r ⊥ r ⊥ ∆ 3 K 1 (∆ 1/2 where the (anti)quark-to-photon fragmentation function at the scale µ F is defined by Eq. (72). This can equivalently be written in terms of the sum of the squared lightcone wavefunctions for the two transverse polarization states of the virtual photon as where the wavefunctions are defined as [41] Ψ T α β (q − , z q , r ⊥ ) = 2π The similarity with the results obtained in Eqs. (76) and (77) for the longitudinally polarized virtual photon can be clearly seen. We will now discuss the computation of the direct photon+quark production cross-section. 24 It is easy to check that the term proportional to im jn J ij T,3 J mn * T,3 is zero because of the symmetric nature of J ij T,3 J mn * T,3 with respect to the indices i and m. and the three terms can be expressed compactly as H qγ T ;reg−reg = C im;jn qq 2 p,r=1 K ij T,p (k ⊥ + k γ ⊥ ) − K ij T,p (k ⊥ + k γ ⊥ − l ⊥ ) K mn T,r (k ⊥ + k γ ⊥ ) − K mn T,r (k ⊥ + k γ ⊥ − l ⊥ ) + C im;jn qq 2 p,r=1 K ij T,p (k ⊥ + k γ ⊥ ) − K ij T,p (k ⊥ + k γ ⊥ − l ⊥ ) K mn T,r+2 (k ⊥ ) − K mn T,r+2 (k ⊥ − l ⊥ ) H qγ T ;reg−ins = 2κ qq and H qγ T ;ins−ins = σ qq K T ins,2 (k ⊥ + k γ ⊥ ) − K T ins,2 (k ⊥ + k γ ⊥ − l ⊥ )