Photoproduction of three jets in the CGC: gluon TMDs and dilute limit

We study the process $\gamma A \to q \bar{q} g+X$ in the Color Glass Condensate (CGC) effective theory. After obtaining the cross section, we consider two kinematic limits which are encompassed in our result. In the so-called correlation limit, the vector sum of the transverse momenta of the three outgoing particles is small with respect to the individual transverse momenta; the cross section then simplifies considerably and can be written in a factorized form, sensitive to both the unpolarized and linearly-polarized Weizs\"acker-Williams transverse momentum dependent gluon distribution function (gluon TMD). The second limit of the CGC cross section that we consider is the dilute limit, which we obtain after performing a weak-field expansion; we recover a typical linear-regime expression, involving a single unintegrated gluon distribution function. Using numerical simulations of the small-$x$ QCD evolution of the TMDs, we investigate the rapidity dependence of the cross section in the correlation limit.


Introduction
Quantum Chromodynamics (QCD) at high energies has been the subject of intense study since many years. In particular, in certain processes such as deep-inelastic lepton-proton scattering, it is possible to probe the constituents of a proton or nucleus at small values of x, where x is the fraction of the proton's longitudinal momentum carried by the struck parton. Since the basic laws of QCD favor the emission of gluons with arbitrary small energies, at low x the proton structure is dominated by the gluon distribution. The power-like rise of this gluon distribution with decreasing x, as predicted by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution equations [1], is expected to be damped by gluon recombination effects once a sufficiently dense regime is reached, characterized by the dynamically generated saturation scale Q s . The ensuing non-linear low-x evolution equations were established both from the point of view of the dilute projectile [2,3], and from that of the highly dense proton or nucleus [4]. In the latter case, the evolution can be obtained from an effective theory known as the Color Glass Condensate (CGC) [5]. Both approaches are equivalent and are often denoted BK-JIMWLK, after their main authors Balitsky, Kovchegov, Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov and Kovner.
Moreover, in the low-density regime, hence at momenta above the saturation scale, the results from linear BFKL evolution are recovered. In this work, we will use the most general semiclassical description of small-x physics, and refer to it as the CGC framework.
More recently, in a series of papers [6] the CGC was explored in reactions characterized by two strongly ordered scales, let us call them P and q: √ s P q ∼ Q s Λ QCD . (1.1) The scales are required to be much smaller than the very high collision energy √ s, but (a priori) larger than Λ QCD and thus perturbative. A prime example for such a reaction is the forward production of a parton pair in lepton-proton or proton-proton collisions, in the so-called correlation limit where the total transverse momentum of the produced pair q ∼ |k 1 + k 2 | is much smaller than the typical momentum of the individual particles P ∼ |k 1 | ∼ |k 2 |. The crucial observation that was made in Refs. [6] is that in these kinematics, the small-x limit of the proton's or nucleus' transverse momentum dependent gluon distributions (gluon TMDs) is probed [7]. Remarkably, at least for 2 → 2 processes, taking the low-x limit of the calculation in the TMD framework yields the same result as taking the correlation limit of the CGC calculation.
For the processes that have been studied, the CGC therefore generalizes two different QCD regimes at low x: the TMD region with two ordered hard scales on the one hand, and the BFKL regime at low enough gluon densities on the other. Inspired on this, a computation scheme dubbed 'Improved TMD factorization' (ITMD) has been developed in Refs. [8], which is applicable to massless 2 → 2 processes and just like the CGC interpolates between the low-x TMD and BFKL (sometimes called High-Energy Factorization, HEF) regimes. Its advantage, however, is that it is much more suitable for numerical implementation. The proof of this scheme as an all-order kinematic twist resummation was recently established by some of us [9].
In addition to the conceptual interest, understanding the TMD framework within the CGC makes it possible to apply the CGC machinery to TMDs, and vice versa. For instance, in [10][11][12] JIMWLK evolution was applied to unpolarized and linearly polarized gluon TMDs, using the nonperturbative McLerran-Venugopalan model [13] as an initial condition. Likewise, Sudakov logarithms, which govern the TMD evolution [14], can be resummed in a consistent way at low x, as was shown in [15] and carried out in e.g. [16].
At least to leading-order accuracy, the above mentioned correspondence between the CGC and TMD frameworks in their overlapping region of validity holds for 2 → 2 processes, also when masses are included [10,12,17], see [18] for an overview. The present work is part of our efforts to investigate whether this is also true for more complicated processes 1 , and can be regarded as a follow up of Ref. [19], in which we calculated the correlation limit of three final state particles (two jets and a photon) in proton-nucleus collisions 2 . In this configuration, the total transverse momentum of the outgoing particles q T = k 1 + k 2 + k 3 is again required to be much smaller than their individual transverse momenta 1 We should remark that recently there was also made progress in this direction from the point of view of TMDs, see Ref. [20]. 2 The correlation limit of two jets and a collinear photon was obtained in [21].
(k 1 , k 1 , k 3 ). However, in contrast to the production of two final state particles, the 2 → 3 kinematics allow us to identify not one but two small transverse sizes in coordinate space. We showed that it is still possible to take the correlation limit, although the procedure becomes more involved, and that once again one arrives at an expression which is factorized in terms of TMDs. We will show in this paper that the same argument holds for three-jet photoproduction, with the additional feature that, due to the simple color structure of the process (since there is no initial-state radiation), only the unpolarized and linearly polarized gluon TMD of the Weizsäcker-Williams type play a role. Using the results from [12], where the JIMWLK evolution of these TMDs was performed, allows us to do a numerical study. The remainder of the paper is organized as follows. In section 2, we give a concise derivation of the CGC cross section for the process γA → qqg + X, relegating most of the details to the appendices. After that, we analytically calculate the correlation limit of our CGC result in section 3, and numerically study its rapidity evolution in section 4. In section 5, we expand our CGC cross section in the dilute limit to recover the high-energy factorization result, and finally summarize our conclusions and outlook in section 6.

CGC cross section
The usual approach to calculate a forward particle production cross section in the CGC is to go to a dipole frame, which is justified at large enough energies and in which the perturbative 'dressing' of the photon state takes place long before the scattering with the highly boosted proton or nucleus target. In this frame, due to Lorentz contraction the projectile sees the target as a highly dense shockwave, and the scattering between both takes place almost instantaneous at a light-cone time which, without losing generality, we can set to be x + = 0. Due to the high energy, it is then justified to describe the scattering in the eikonal approximation, where the hard particles do not undergo a change in transverse position upon interacting with the target's gluon fields, but only exchange transverse momentum and color. Moreover, since due to the high density these gluon fields are semiclassical, i.e. g s A ∼ 1, the interactions need to be resummed which is done by the use of Wilson lines. Finally, the way one averages over the semiclassical target gluon fields reflects the nonperturbative proton structure, which we leave unspecified.
In light-cone perturbation theory (LCPT), the partonic cross section is defined as the expectation value of the number operator calculated in the relevant component of the outgoing photon wave function: where the factor 1/2 on the right hand side stems from averaging over the photon polarization λ. The vectors k i ≡ (k + i , k i ) stand for the three-momenta of the produced particles, with k + i being the longitudinal and k i the transverse momenta ( k 1 for the quark, k 2 for the gluon and k 3 for the antiquark), while p + is the forward component of the incoming photon (note that we assume p = 0). Moreover, the number operators N i for the dressed (this will be important later on) quark, gluon and antiquark, are defined in terms of the corresponding creation and annihilation operators. For example, the gluon number operator N g is defined as ) the creation (annihilation) operator of a dressed gluon with color a, polarization i and three-momentum k 2 . The action of these operators on dressed Fock states is defined in the standard way. For example, for the gluon we have In order to compute the partonic cross section, we need to calculate the explicit expression of the outgoing Fock state to the relevant order. Hence, for the process in this work, one needs the outgoing wave function for a quark-antiquark-gluon final state initiated by a real photon. Since the derivations of similar outgoing wave functions, i.e. q → qgγ and g → qqγ, have been presented in detail Ref. [19], we do not show its derivation here but rather give a summary in the appendices A.1 and A.2. The order g s g e contribution to the dressed outgoing wave function of a real photon with longitudinal momentum p + and vanishing transverse momentum is given by 3 : In the above, the quark with color i and spin s carries longitudinal momentum k + 1 and is located at the transverse position x 1 . Likewise, the gluon with color c and polarization η has a longitudinal momentum k + 2 and is located at the transverse position x 2 , and finally the antiquark with color j and spin s carries longitudinal momentum k + 3 = p + − k + 1 − k + 2 and sits at the transverse position x 3 . The longitudinal momentum fractions that enter the splitting functions are defined as The outgoing state Eq. (2.4) contains combinations of three splitting functions, correspond- 3 We have introduced a shorthand notation for the two dimensional coordinate integrals: x = d 2 x. ing respectively to the splitting of a photon into a quark-antiquark pair, the radiation of a gluon by a quark, and the instantaneous splitting of a photon into a gluon, quark and antiquark. Their explicit expressions are: where σ 3 is the third Pauli matrix, and with ij the Levi-Civita tensor in two dimensions, with 12 = +1. Note that the splitting function for an antiquark emitting a gluon, is obtained from the corresponding splitting function for the quark by an overall multiplication with −1 and exchanging the spin indices. The final ingredients of Eq. (2.4) are the functions G q , Gq and G C , which incorporate the Wilson line structures of the three different amplitudes where the gluon is emitted from the quark, from the antiquark, or directly from the photon during its instantaneous splitting into a quark-antiquark-gluon state (see Fig. 2.1). They are defined as: In these functions, the Wilson lines are defined in the standard way as the path ordered exponential of the semiclassical gluon field α − a (x + , x) of the target in a covariant gauge: where t a is the SU (N c ) generator in the fundamental representation, indicated with the subscript 4 F . Moreover, in the above expressions, A i (x) stands for the standard non-Abelian Weizsäcker-Williams field, which accounts for a single emission: On the other hand, A i (ξ, x; χ, y) is the modified Weizsäcker-Williams field, responsible for two successive emissions, whose explicit expression is written as (2.12) Finally, we have introduced the Coulomb field which accounts for the instantaneous emission, defined as: (2.13) As mentioned earlier, the partonic cross section Eq. (2.1) is formally defined as the expectation value of the number operator in the outgoing wave function. This expectation value can be calculated in a straightforward way from the action of the creation/annihilation operators on the outgoing wave function given in Eq. (2.4). The resulting cross section is then averaged over the configurations of the target field α − (x + , x). We denote this averaging procedure by · · · x A , since it introduces an implicit dependence on x A : the longitudinal momentum fraction of the gluons in the target wave function. The result is organized in the following way: (2.14) Each contribution in Eq. (2.14) can be studied separately. Let us start with what we call the 'quark-quark' contribution I qq , corresponding to a gluon emission from the quark (leftmost panel in Fig. 2.1) both in the amplitude and in the complex conjugate amplitude: where the function Mλλ ;ηη qq is the product of the γ → qq and q → qg splitting functions, defined in Eq. (2.6), and reads It can be evaluated explicitly by using the expressions in Eq. (2.6), yielding: The color trace in Eq. (2.15) can be performed by using the explicit expression of the function G q , given in Eq. (2.7), and using the Fierz identity The result of the color algebra can be written in a convenient way by introducing the following functions where the dipole and the quadrupole operators are defined as The final expression for the I qq contribution to the partonic cross section is: The second contribution to the partonic cross section Eq. (2.14) is the 'antiquarkantiquark' term Iqq, corresponding to the emission of the gluon from the antiquark (see the middle panel in Fig. 2.1) both in the amplitude and in the complex conjugate amplitude. It reads: In order to arrive to the explicit expression of the antiquark-antiquark contribution, one computes the trace over color indices and performs the color algebra by using the explicit expression of the function Gq given in Eq. (2.8 Similarly to the quark-quark contribution, the function Mλλ ;ηη qq encodes the products of the γ → qg andq →qg splitting functions: Note that, due to the charge conjugation symmetry on the cross section level, Mλλ ;ηη qq can be obtained from M λλ ;ηη qq by swapping the quantum numbers of the quark and the antiquark. Alternatively, we can calculate it using the explicit expressions of the splitting functions given in Eq. (2.6), yielding: The next contribution that we consider is I qq , corresponding to the emission of the gluon from the antiquark in the amplitude and from the quark in the complex conjugate amplitude, or vice versa. This term can be written as: (2.28) As in the previous two cases, to evaluate this contribution we use the explicit expressions for the functions G q and Gq, given in Eqs. (2.7) and (2.8), respectively. This yields: where we have introduced a new function W 4 x 1 , x 2 ; y 1 , y 2 corresponding to the following combination of dipole operators: The product of splitting amplitudes Mλλ ;ηη qq for this contribution is: and reads, when computed explicitly: The next step is to compute the contributions stemming from the instantaneous splitting of the photon into a quark, an antiquark and a gluon (rightmost panel in Fig. 2.1). Let us start with the I CC contribution, which corresponds to the instantaneous splitting both in the amplitude and in the complex conjugate amplitude. This contribution is written in terms of the function G C , defined in Eq. (2.9), and reads After performing the color algebra, one obtains: The product of splitting functions for the I CC contribution reads The remaining two contributions that need to be computed are I Cq x A and I Cq x A , corresponding the interference between the gluon emission from the quark resp. antiquark in the amplitude, and the instantaneous splitting in the complex conjugate amplitude. The first one is given by: As in the case of other contributions, using the explicit expressions of the functions G q and G C and performing the color algebra results in the following expression: with the product of splitting functions Mλη Cq equal to Finally, the last contribution I Cq x A is written in terms of the functions G C and Gq as follows: and yields, after performing the color algebra: where the product of splitting functions Mλη Cq is: Let us summarize our findings for this section. The partonic cross section for the photoproduction of a quark, an antiquark, and a gluon, is given in Eq. (2.14). Each contribution to this cross section is calculated separately and the final results are given in Eqs. (2.23), (2.25), (2.29), (2.34), (2.37) and (2.40). As we would expect from the charge conjugation symmetry of QCD, the result is fully symmetric under the exchange of the quark with the antiquark.
Finally, the photoproduction cross section in ep or in eA collisions, or alternatively in ultra-peripheral pp or pA collisions (UPCs), can be obtained in a very straightforward way from the partonic cross section by using the equivalent photon approximation. This approximation consists in simply convolving the partonic cross section with the relevant photon flux f e,P,A→γ (y, µ 2 ): where y = P · p/P · p + / + is the longitudinal momentum fraction carried by the photon (with µ the momentum of the photon source) and where µ 2 is the factorization scale. For example, if the photon source is an electron, the real photon flux f e→γ (y, µ 2 ) is given in the well-known Weizsäcker-Williams approximation by the formula [22]: where m e is the electron mass and α the fine-structure constant.

Correlation limit and gluon TMDs
In recent years, there has been a lot of activity in the study of gluon TMDs from CGC calculations. In Refs. [6], it was shown that the dijet production cross section in forward pA collisions can be written in terms of gluon TMDs, in a specific kinematic limit that is referred to as the correlation limit. In these specific kinematics, the total transverse momentum of the produced pair (∼ |k 1 + k 2 |) is assumed to be much smaller than the typical momentum of the jets (∼ |k 1 | ∼ |k 2 |). In such a situation, the produced jets fly almost back-to-back in momentum space, which in coordinate space corresponds to a small transverse distance between the produced jet pair, usually referred to as the dipole size.
Recently, the correlation limit of three final state particles (two jets and a photon) has been studied in Ref. [19]. In this configuration the total transverse momentum q T = k 1 + k 2 + k 3 is again required to be much smaller than the individual transverse momenta of the produced jets (k 1 , k 1 , k 3 ). In contrast to the production of two final state particles, these kinematics allow us to identify not one but two small transverse sizes in coordinate space.
Let us turn to our process, and identify the small transverse sizes around which we can perform a Taylor expansion. These sizes are not necessarily the same for each subprocess, and can be identified by inspecting the denominators of the structures Eqs. (2.11), (2.12) and (2.13), which appear in the functions G q , Gq and G C defined in Eqs. (2.7), (2.8) and (2.9), respectively.
Firstly, let us consider the function G q , Eq. (2.7), which is a part of the amplitude where the final state gluon is emitted from the quark. In the correlation limit, the small transverse distance for this contribution appears to be where r g is identified as the size of the dipole formed by the final quark and gluon, and where rq is the dipole size of the final state antiquark and the intermediate quark.
With the help of the delta function δ (2) x 2 , which accompanies the function G q in the outgoing photon wave function Eq. (2.4), we can substitute the coordinates x 1 , x 3 and v in favor of the above defined dipole sizes: The small-dipole expansion of G q is then straightforward, and reads: . Likewise, the dipole sizes appearing in the function Gq, Eq. (2.8), which corresponds to the emission of the final state gluon from the antiquark, are: where r q is the dipole size of the final quark and the intermediate antiquark, and where in contrast to the previous case r g is dipole size of the final gluon and anti quark. Taking the delta function δ (2) x 3 into account, which accompanies the function Gq (see Eq. (2.4)), the transverse coordinates can be written as: The Taylor expansion of Gq then yields: Finally, the small transverse sizes that appear in the function G C , Eq. (2.9), corresponding to the instantaneous emission of the quark-gluon-antiquark from the incoming photon, are: with r g the transverse size of the final quark-gluon dipole, and rq the size of the dipole formed by the final state antiquark with the incoming photon. Again, rewriting the transverse space coordinates in terms of these dipole sizes we get: Then, the Taylor expansion of the function G C reads With the small-dipole expansions of the functions G q , Gq and G C at hand, we are now ready to compute the correlation limit of the partonic cross section Eq. (2.14). Let us start with the I qq contribution, given in Eq. (2.15). Using the change of variables that are introduced in Eq. (3.2) both in the amplitude and in the complex conjugate amplitude, as well as the Taylor-expanded expression of the function G q given in Eq. (3.3), we obtain the following result: where the combinations of transverse momenta q T and Q are defined as 11) and are ordered as |q T | |k 3 |, |Q|, corresponding to the dipole sizes to which they are conjugate being ordered as |x 2 |, |x 2 | |r g |, |r g |, |r q |, |rq|. The integrations over the dipole sizes r g , r g , rq and r q factorize from the Wilson line structure, and can be performed with the help of the following two integrals 5 : and rgrq e −iQ·rg−ik 3 ·rq Aη(r g ) r ī q Aλ ξ 3 , rq; where we introduced (3.14) Using the results Eqs. (3.12) and (3.13), Eq. (3.10) can be simplified to: where the hard factor is defined as and where we introduced the compact notation: The crucial observation is now that the remaining integrals of x 2 and x 2 in Eq. (3.15), over the Wilson line structure, are nothing but the small-x limit of the so-called Weizsäcker-Williams (WW) gluon TMDs 6 (see e.g. Ref. [12]): The derivations of these two integrals can be found in [19]. 6 Note that in our earlier work, these TMDs are usually written as F gg and H gg .
In the above, F gg (x A , q T ) is the unpolarized Weizsäcker-Williams gluon TMD, and H gg (x A , q T ) its linearly polarized partner. Substituting Eq. (3.18) into Eq. (3.15), we can write the final result for I qq x A in the following TMD-factorized form: Let us now calculate the correlation limit of the Iqq x A contribution to the partonic cross section, which corresponds to the gluon emission from the antiquark both in the amplitude and in the complex conjugate amplitude. The calculation is straightforward and can be performed following the same steps as in the previous case. Equivalently however, we can circumvent the calculation by observing that, since QCD preserves Cparity, Iqq x A has be equal to I qq x A after exchanging the quark with the antiquark. In our notation, this corresponds to exchanging 1 ↔ 3, as well as swapping the relevant color indices (note that, with the introduction of the products of splitting functions M, the spin indices are already contracted). Hence, we have: where the new combination of transverse momentum K is defined as and |q T | |k 1 | ∼ |K|. Thus, one can immediately write down the following result: where the hard part is defined as The interference contribution I qq x A originates from the gluon emission from the antiquark in the amplitude and from the quark in the complex conjugate amplitude, or vice versa. The small transverse dipole sizes that appear in this contribution are given by Eqs. (3.2) and (3.4). When written in terms of these new variables, this contribution takes the following form: After introducing the expanded expressions for the functions G q and Gq given in Eqs. (3.3) and (3.6), respectively, one can easily perform the transverse integrals. The final expression for the I qq x A contribution in the correlation limit reads: where the hard part is equal to: We can now continue with those contributions to the cross section that include the instantaneous emission of the quark-antiquark-gluon final state from the incoming photon. The first of these contributions is I CC x A , corresponding to the instantaneous emission both in the amplitude and in the complex conjugate amplitude. After introducing the change of variables given in Eq. (3.7), and using the expanded expression of the function G C , Eq. (3.9), I CC x A can be written as: The integrals over the dipole sizes r g , r g , rq and r q can be performed thanks to the following expressions Finally, one obtains: with the hard part equal to The last two contributions to the cross section are the terms with the instantaneous emission in the conjugate amplitude, and the gluon emission from the quark resp. antiquark in the amplitude. Their calculation is completely analogous to the previous cases, and for I Cq x A one obtains: Likewise, the result for I Cq x A is: where the hard part reads: Combining all separate contributions, we reach the main result of this section, i.e. the TMD factorized expression for the photoproduction of three jets in the correlation limit at low-x: where F W W (x A , q T ) and H W W (x A , q T ) are the unpolarized and linearly polarized Weizsäcker-Williams gluon TMDs defined in Eq. (3.18), and where the hard factor is found to be:

Numerical study of the cross section in the correlation limit
In this section, we will further study the cross section in the correlation limit, using both analytical and numerical models of the gluon TMDs in the target proton. Before doing so, let us make some estimates of the relevant phase space for our process. Photoproduction can take place as the low-Q 2 limit of deep-inelastic lepton-proton/ion scattering, for instance in the future Electron-Ion Collider, or in ultra-peripheral collisions (UPCs) involving protons or heavy ions. Let us denote our process as l( ) + A(p A ) → q(k 1 ) + g(k 2 ) +q(k 3 ) + X, with A the target proton or ion, and where l the source of the real photon flux; either a charged lepton, either a proton or ion. The momenta of the photon source and the target are: and hence the center-of-mass energy is equal to s = ( + p A ) 2 2 + p − A . The four-momenta of the real photon p and the gluon k are then given by and the value of x reached in the experiment is: It follows that for realistic EIC center-of-mass (c.o.m.) energies √ s ∼ 100 GeV, and for values of y as large as possible, demanding that x 10 −2 means that the transverse momenta of the jets should not exceed 3 GeV. On the other hand, in ultra-peripheral proton-proton or lead-lead collisions at the LHC one is able to reach c.o.m. energies of the order of √ s ∼ 3 TeV and √ s ∼ 7 · 10 2 GeV, respectively (see e.g. Ref. [26]). In such collisions, values of x down to x ∼ 10 −3 − 10 −4 are attainable, with transverse jet momenta of more workable magnitudes |k i | ∼ 10 GeV.
For these kinematics, it is a reasonable choice to use the nonperturbative McLerran-Venugopalan (MV) model [13] to estimate the gluon distributions inside the unpolarized nucleus or proton. In the MV model, the unpolarized and linearly polarized Weizsäcker-Williams gluon TMDs read [6]: , .

(4.4)
In the above formulas, S ⊥ is the transverse size of the proton, and Λ is an infrared cutoff, which we take to be Λ QCD 0.2 GeV). Furthermore, Q sg is the gluon saturation scale for which we choose the value Q sg = 0.6 GeV for x = 10 −2 [11,12]. A factor e is added inside the logarithms to guarantee the convergence of the expressions Eq. (4.4).
With these expressions at hand, we plot the logarithm of the differential real photonproton cross section in the correlation limit, Eq. (3.35), but expressed in slightly different variables and with the overall azimuthal dependence integrated out: ln dσ γA→qqg+X dθ 12 dθ 32 Π i dξ i d|k i | corr. limit . (4.5) In the above, the angles θ 12 and θ 32 are defined as in Fig. 4.1. We present the result in Fig. 4.2 for the choice |k i | = 10 GeV and ξ i = 1/3 with i = 1, 2, 3, and using the WW gluon TMDs introduced in Eq. (4.4). Since our calculation is valid in the correlation limit |q T | |k i | ∼ 10 GeV, the maximum value of q T where our calculation can be trusted is estimated to be q T ,max = √ 10 GeV, and is illustrated by a black dashed line. Note that we choose to normalize the cross section by its maximal value over the kinematic range, so that we are not sensitive to prefactors such as S ⊥ . For illustrative purposes, we separately plot the contributions to the cross section from the unpolarized and from the linearly polarized gluon TMDs. Since the amplitude of the latter is several orders of magnitude smaller than the former, its contribution to the total cross section is almost negligible. Nevertheless, as is clearly visible from the plots, the contribution from the linearly polarized gluon TMD exhibits azimuthal modulations, which can be exploited to extract this TMD experimentally (see e.g. Ref. [23]).
Furthermore, in Figs. 4.3 and 4.4, we show the nonlinear low-x evolution of the cross section. For this, we make use of the results of Ref. [12], in which the numerical JIMWLK evolution of different gluon TMDs was performed, including the unpolarized and linearly  polarized Weizsäcker-Williams distribution. Its initial condition is a numerical implementation of the MV model, which slightly differs from the analytical one in Eqs. (4.4). Once again, we identify the starting point of the evolution to be x = x 0 = 10 −2 with an associated saturation scale Q sg = 0.6 GeV. The evolution is performed to values of x 10 −3 and x 10 −4 , with corresponding saturation scales of Q sg = 0.86 GeV and Q sg = 1.5 GeV, respectively [11,12]. Clearly, the peak of the cross section around q T = 0 broadens quickly with the evolution, which is expected from the behavior of the saturation scale Q sg around which the gluons inside the target are distributed. Likewise, the angular modulations, which appear in the contribution to the cross section from linearly polarized gluons, are suppressed and pushed further away from the center at q T = 0.

Weak-field approximation in the CGC and the HEF limit
As is well known, the high-energy factorization (HEF) framework is compassed within the CGC formalism, and can be extracted by performing the weak-field approximation for the gauge fields of the target. In this limit, the Wilson line structure in the amplitude and in the complex conjugate amplitude should be expanded to second order in powers of the background field α − a (x + , x) of the target. However, the structure of the functions G q , Gq and G C , which incorporate the Wilson lines at the amplitude level, ensure that it is enough to perform the expansion to the first order in powers of the background field since the zeroth order terms vanish (see Eqs. (2.7), (2.8) and (2.9)). The standard weak-field expansion of a fundamental Wilson line reads In order to compute the dilute limit of the cross section, we will follow the same strategy as in the correlation limit, i.e.: we first perform the expansion on the level of the amplitude for the quark, antiquark and instantaneous contributions, and then combine the results on the cross section level. To do so, we introduce the reduced quark amplitude As is clear from the above expression, we define the reduced amplitude as the regular amplitude, stripped from the tensor structure of the product of the splitting functions Mλλ ;ηη , such that, taking the example of the I qq contribution to the cross section (2.15): After introducing the following change of variables and using Eq. (5.1) for the expansion of the Wilson lines in the function G q , the reduced quark amplitude can be cast into the following compact form: In the above expression, q T is the total transverse momentum defined in Eq. (3.11), and we defined the following combinations of transverse momenta: In Eq. (5.5) the integrals over r andx 3 are factorized from the rest of the expression and can be performed in straightforward manner. The reduced quark amplitude then becomes: whereÃ λ (k) andÃ λ (k; c, p) are the standard and modified Weizsäcker-Williams fields in momentum space, given by: (5.8) and the combination of forward momentum fractions c 0 = (ξ 1 ξ 2 )/(ξ 2 3 ξ 3 ) which was defined earlier.
In a similar way, the reduced antiquark amplitude is defined as Following the same steps as for the reduced quark amplitude, the weak-field limit of the reduced antiquark amplitude can be computed and reads: 11) and the combinationc 0 = (ξ 2 ξ 3 )/(ξ 2 1 ξ 1 ) which was also defined earlier. Finally, we define the reduced amplitude for the instantaneous contribution: which, after the weak-field expansion, takes the following form where the different combinations of transverse momenta p 13 , p 23 and m 12 are defined in Eq. (5.6). The functionC(ξ 3 ,ξ 3 p; ξ 1 ξ 2 , k) is the Fourier transform of the Coulomb field defined in Eq. (2.13) and reads: (5.14) We will now use the explicit expressions of the reduced amplitudes in the weak-field limit: Eq. (5.7), (5.10), and (5.13), to calculate the HEF expression for the photoproduction cross section.
Let us start with the I qq x A contribution to the partonic cross section. Written in terms of the weak-field expansion of the reduced quark amplitude Eq. (5.7), this contribution can be cast into the following form: where the Weizsäcker-Williams field structure Wλλ ;ηη qq for this contribution is computed as Using the relation (A.72) between the correlator of gauge fields and the unintegrated gluon PDF F g/A , we finally obtain: The Iqq x A contribution can be computed in a similar manner. Using the dilute expansion of the reduced antiquark amplitude given in Eq. (5.10), together with the definition of the unintegrated gluon distribution, one obtains: Not surprisingly, the above expression is the same as the one we obtained for the quarkquark contribution in Eq. (5.16), up to the change of the subscripts 1 ↔ 3. Since, as discussed previously, the functions that encode the product of splitting amplitudes follow the same rule, we have that 20) reflecting the invariance of QCD under charge conjugation. It turns out that the HEF limits of all the remaining contributions, i.e. I qq x A , I CC x A , I Cq x A and I Cq x A , which can be calculated with the help of the expansions Eqs. (5.7), (5.7) and (5.13) on the amplitude level, can be cast in a similar form: where the symbol ⊗ denotes a contraction over all open vector indices. The corresponding Weizsäcker-Williams field structures W are given by, respectively: and lastly: We can now combine all the contributions that have been computed separately, and write down the final factorized expression for the photoproduction of three jets in the HEF limit as: where the total Weizsäcker-Williams field structure W total reads W total = π 2 N c g 2 e g 4

Summary and outlook
We calculated the partonic cross section for the forward inclusive production of a quarkantiquark pair plus a gluon in the scattering of a real photon with a proton or nucleus. The computation was performed at leading order in the CGC effective theory, valid at low values of x, taking the multiple rescatterings of the partons off the semiclassical gluon fields in the target proton or nucleus into account. From our partonic cross section, the cross section for three-jet production in the low-Q 2 limit of deep-inelastic scattering, or in ultraperipheral collisions, can be easily obtained by convolving the result with the relevant real photon flux. Note that the cross section for three-jet production in deep-inelastic scattering at low-x has been calculated in a helicity framework in [24].
Our main result is the correlation limit of this cross section, corresponding to the kinematical configuration where the total transverse momentum q T of the three outgoing particles is much smaller than their individual transverse momenta. We demonstrated how in this limit the cross section simplifies and factorizes into a partonic hard part on the one hand, independent of q T , and a gluon correlator on the other hand, which parametrizes the proton or nucleus target in terms of the unpolarized and linearly polarized gluon TMDs F W W (x A , q T ) and H W W (x A , q T ). Using earlier results of the JIMWLK evolution of these two gluon TMDs, we numerically studied the nonlinear small-x evolution of this cross section.
In addition, we calculated the dilute limit of the CGC cross section by performing the weak-field expansion. Once again, a simple factorization formula was obtained, this time in terms of the unintegrated gluon PDF, and with a hard part corresponding to the γg * → qqg off-shell matrix element in the language of high-energy factorization (or k T -factorization).
Our calculation provides an important contribution towards establishing the phenomenology of forward particle production at low-x, in the line of [20,24,32]. In particular, recently there has been much interest in processes with photons in the final state as a clean probe for saturation effects [19,[33][34][35][36], and it can be noted that the cross section for the photoproduction of two jets and a photon can very easily be adapted from the results of the present paper, with a simplified Wilson line structure.
More importantly, we demonstrated for the first time that the CGC-TMD correspondence still holds, at least at leading order, when considering three outgoing colored particles. Finally, within the framework set up in this work and in [19], we are currently pursuing the computation of photon-jet production in the CGC at NLO, to investigate whether this correspondence holds beyond tree level as well.
In the above expression, the standard properties of light-cone perturbation theory (LCPT) apply. In particular, at each vertex three-momentum is conserved, hence it is understood that q 1 = p − q 2 , k 3 = p − k 1 − k 2 , and so on. Moreover, the momenta satisfy on-shell conditions, implying that in our massless case Introducing what we call the 'wave functions' F , the dressed photon state can be cast in a much more compact form: The wave functions contain the dynamics of the splittings, and their expressions in terms of the matrix elements are obtained by comparing Eqs. (A.2) and (A.3). They can be calculated with the help of the Feynman rules in appendix A.3. Since their calculation has been performed explicitly for similar processes in Refs. [19,21], but with slightly different conventions, we restrict ourselves here to the computation of F (1) γ : the leadingorder splitting of a photon into a quark-antiquark pair. For the other wave functions, we will merely present the final results.
The (dimension −1) γ → qq wave function (see Fig. A.1) is defined as: where one should be careful to note that an arrow over a momentum vector indicates k = (k + , k), while an arrow over a coordinate means x = (x − , x). Making use of the LCPT Feynman rules as well as the conventions for the (anti-)commutation relations listed in appendix A.3: Note that, with a minor abuse of notation, we extracted the factor g e δ ij from the interaction Hamiltonian and placed it in front of the wave function in Eq. (A.3). After some algebra, we find that: such that we obtain: As we will show in appendix A.4, the numerator in the above formula can be calculated in terms of the so-called good spinors, on which we will elaborate later. In the massless case under consideration, the result is generic and does not depend on whether one considers quark or antiquark spinors. We can therefore read off the result from the general formula (A.37), and adapt it to the kinematics of the γ → qq splitting under consideration: To make this result more explicit, we choose a set of polarization vectors and Dirac spinors, for which we refer to appendix A.6. With this choice, our final result for the γ → qq wave function is: where Ψ λλ ss (ξ) is the dimensionless splitting function, defined as: For the other wave functions, the computation is analogous, and we refer to Refs. [19,21] for the explicit calculation of similar second-order functions. Apart from the γ → qq case, the q → qg andq →qg leading-order wave functions (see Fig. A.1) will be important in this work: Once again, with some abuse of notation the factors g s t c ij were extracted. The results for the second-order splitting functions (see Fig. 2.1) are: where the notations ξ 1 = k + 1 /p + , ξ 2 = k + 2 /p + , ξ 3 = (p + − k + 1 − k + 2 )/p + were introduced. The splitting function of the instantaneous splitting of a photon into a gluon, a quark and an antiquark, is given by: In mixed Fourier space, the wave functions are found to be: and The expressions for the non-Abelian Weizsäcker-Williams field A, the modified Weizsäcker-Williams field A and the Coulomb field C can be found in Eqs. (2.11), (2.12) and (2.13).

A.2 Outgoing Fock state
The first step towards the computation of the outgoing Fock state, is performing the Fourier transform of the transverse coordinates. For the zeroth-order term of the dressed photon state (A.3), this is trivial and yields: of SU (N c ) are in the fundamental or adjoint representation, respectively): Finally, the outgoing state needs to be written as a function of the dressed states, not the bare ones. Schematically, the procedure goes as follows: first, observe that all the dressed states are, up to g e g s accuracy, related to the bare ones as: The outgoing state is then given by: The above procedure can be performed explicitly from a closer look at the |qq → |qqg splittings. It is easy to see that, in mixed Fourier space: We thus obtain for the outgoing state (retaining only the terms that yield the three outgoing particles we are after): one arrives at the final expression Eq. (2.4).

A.3 LCPT conventions and Feynman rules
We follow the conventions of Ref. [25], in which the quark and gluon fields are defined as follows in terms of the creation-and annihilation operators: with the following (anti-)commutation relations: As is clear from the above definitions, in our conventions all the creation and annihilation operators have mass dimension −1. Hence, it follows that the for the moment unspecified Dirac spinors are required to have dimension 1/2, in order for the quark field to have dimension 3/2. Following the same reasoning, the polarization vectors are dimensionless, such that the boson field has dimension 1. With these definitions, each n-particle Fock state in momentum space has dimension −n. As in the covariant case, the quark-antiquark-photon vertex is simply obtained from the quark-antiquark-gluon one by replacing g s → g e and A g → A γ . Note that : O : denotes normal ordering, and that we used the notation x = (x − , x).
The above conventions are in contrast with the ones used in our earlier work Ref. [19,21], where all creation and annihilation operators have dimension −3/2, while the quark and gluon fields keep their usual dimensions 3/2 and 1.
The Feynman rules for LCPT can be found in Refs. [25,28]. The relevant ones for our process are listed in Fig. A.2, where it is understood that the gluon field is A g ≡ A g,c t c . The quark-antiquark-photon vertex is then obtained from the quark-antiquark-gluon one by the replacement g s → g e and A g → A γ .

A.4 Generic quark-boson-antiquark splitting function
In this subsection, we show how to derive a generic expression for the quark-bosonantiquark splitting, irregardless of which particle is the radiator. Without losing generality, let us focus on the structure:ū ( q 2 ) ¡ λ ( p)v( q 1 ) . (A.31) We will decompose the spinors in terms of good and bad components, which can be performed by means of the projectors: In LCPT even intermediate particles follow the classical equations of motion, hence it follows in a straightforward way from the Dirac equation that the bad components of a spinor depend on the good ones through the relations (see Ref. [25]): and similarly:ū With the above relations, as well as the following parameterization of the polarization vector for the gluon: and using the identities in appendix A.5, we obtain the following expression for Eq. (A.31): (A.36) Furthermore, using the decomposition Eq. (A.42), we obtain the result: (A.37) which was obtained earlier in Ref. [25]. Since in the massless case the good and bad projections (A.33) and (A.34) are the same for the quark and the antiquark, the above result is general, and one can freely interchangeū →v or v → u.

A.6 Explicit spinor representation
One can simplify the wave functions considerably by choosing an explicit parameterization of the Dirac spinors. We follow the so-called 'Kogut-Soper' conventions [27,28], in which the transverse polarization vectors of the gluon are written as follows: With these conventions, the Dirac spinors are given by (1 corresponds to spin up, −1 to spin down): and v 1 ( k) = 1 (A.53) From the above expressions and in the massless case, it is very easy to obtain the good spinors u s G ( k) = P G u s ( k) and v s G ( k) = P G v s ( k):

A.7 Weizsäcker-Williams gluon TMD and unintegrated gluon distribution
In this part we construct the gluon distributions of a right-moving nuclear target, to be consistent with the rest of the appendix in which we set up the LCPT rules for a right mover. However, in the calculation to which this paper is devoted, we take the projectile -the dressed photon-to be a right mover, scattering off a left-moving target. In order to apply the formulas derived in this subsection to the calculation in the main text, they should be converted to a left-moving target which can be done by simply exchanging the indices + ↔ −. Note that, in this appendix, we also keep the gluon fields A µ general; i.e. they are not necessarily semiclassical background fields α µ . We define the gluon density inside a right-moving proton or nucleus as the gluon number operator evaluated in the target's nonperturbative Fock state [29]: where x A = k + /p + A is the fraction of the target's longitudinal momentum carried by the gluon. From the Fourier expansion of the gluon quantum field, keeping in mind that the integration over the longitudinal momentum component is restricted to k + > 0, we easily obtain that: where we chose linear polarization vectors i λ = δ i λ . Moreover, in a light-cone gauge A + = 0 the color field strengths are equal to F i+ c ( k) = −ik + A i c ( k), such that: Weizsäcker-Williams gluon TMD The operator definition of the Weizsäcker-Williams gluon TMD is These links can be eliminated by choosing a light-cone gauge 7 A + = 0, A i (ξ − = +∞) = 0, such that, where we also evaluated the trace, using Tr(t a t b ) = δ ab /2. Furthermore, making use of translational invariance and Fourier transforming the field strengths: Finally, the Fock states of the nuclear target differ from the usual hadronic states by a normalization: such that: Unintegrated gluon distribution The unintegrated gluon distribution F g/A (x A , k), on the other hand, is commonly defined via the dipole cross section in the dilute limit of the CGC, to which it is linearly related [30]: To compute the dipole cross section in the dilute limit (again, for a right-moving target), we need to expand the Wilson lines to second order in the gauge fields, this time in a covariant gauge: Note that, to the present accuracy the path ordering does not matter (since Tr(t a t b ) = Tr(t b t a ) = δ ab /2), and there is thus a strict separation between transverse and longitudinal dynamics. With this expansion, the dipole cross section reads: which is equal to Eq. (A.69) when we define: Restoring the factor e −ik + (v − −w − ) 1, and choosing a covariant gauge such that F i+ a (x − , x) = ∂ i A + a (x − , x), we obtain after partial integration: Hence, in the dilute limit at low-x where Eq. (A.69) holds, the unintegrated gluon PDF and the Weizsäcker-Williams gluon TMD are equal up to a factor π [31]. In the appropriate gauge, both can be shown to be equal (up to the same factor π) to the gluon number counting operator evaluated in the CGC average.