Forward dijets in proton-nucleus collisions at next-to-leading order: the real corrections

Using the CGC effective theory together with the hybrid factorisation, we study forward dijet production in proton-nucleus collisions beyond leading order. In this paper, we compute the"real"next-to-leading order (NLO) corrections, i.e. the radiative corrections associated with a three-parton final state, out of which only two are being measured. To that aim, we start by revisiting our previous results for the three-parton cross-section presented in our previous paper. After some reshuffling of terms, we deduce new expressions for these results, which not only look considerably simpler, but are also physically more transparent. We also correct several errors in this process. The real NLO corrections to inclusive dijet production are then obtained by integrating out the kinematics of any of the three final partons. We explicitly work out the interesting limits where the unmeasured parton is either a soft gluon, or the product of a collinear splitting. We find the expected results in both limits: the B-JIMWLK evolution of the leading-order dijet cross-section in the first case (soft gluon) and, respectively, the DGLAP evolution of the initial and final states in the second case (collinear splitting). The"virtual"NLO corrections to dijet production will be presented in a subsequent publication.


Introduction
In this paper we shall study inclusive dijet production in proton-nucleus (pA) collisions at forward rapidities, i.e. at very small angles with respect to the collision axis in the fragmentation region of the proton. This kinematics is interesting in that it gives us access to the small-x part of the gluon distribution in the nuclear target, where high-density phenomena like gluon saturation are expected to be important. For this particular set-up, it is a good approximation to assume that the final statethe two measured jets possibly accompanied by unmeasured particles -is produced via radiation from a single parton from the incoming proton. The role of the scattering is to put on-shell the otherwise virtual quanta from the parton (light-cone) wavefunction and also to give these quanta a non-trivial distribution in energies and momenta. In particular, their final distribution in transverse momenta (or azimuthal angles) should reflect properties of the nuclear gluon distribution, that we are ultimately interested in. To avoid a proliferation of cases, we shall restrict ourselves to the quark channel, that is, we shall consider only the case where the parton from the proton which participates in the scattering is a quark. The addition of the gluon channel is in principle straightforward and will be addressed in a subsequent publication. As implicit in the previous discussion, we shall use the framework of the so-called "hybrid factorisation" [2][3][4], where the leading quark is assumed to be collinear with the proton (and described by the standard quark distribution), whereas its interactions with the nuclear target are described within the Colour Glass Condensate (CGC) effective theory [5,6]. The nucleus is viewed as a source of strong, random, colour fields representing the small-x gluons and their correlations. The high-energy scattering between a parton from the proton -the leading quark and its radiation -and these strong colour fields will be computed in the eikonal approximation, that is, by associating a Wilson line (describing colour precession) to each parton projectile. As a result, the cross-section corresponding to a given (partonic) final state can be related to gauge-invariant products of Wilson lines, whose high-energy evolution (via soft gluon emissions) is described by the Balitsky-JIMWLK equations [7][8][9][10][11][12][13][14] -actually, an infinite hierarchy of coupled equations for multi-point correlations. This evolution becomes considerably simpler in the limit of a large number N c 1 of colours, where the hierarchy acquires a "triangular" structure. In particular, the first equation in this hierarchy reduces to a closed, non-linear, equation -the Balitsky-Kovchegov equation [7,8] -for the elastic S-matrix of a colour dipole.
The formalism as a whole -meaning the non-linear evolution equations with increasing energy and the hybrid factorisation -has been originally proposed to leading order (LO) in perturbative QCD [4][5][6]. But for realistic applications to the phenomenology, one needs next-to-leading order (NLO) accuracy, at least. Note that, unlike for the collinear factorisation (and the associated DGLAP evolution), the compatibility between the high-energy factorisation and the weak coupling expansion is a a priori unclear, at both conceptual and practical level. Yet, via explicit calculations and tenuous efforts, one was able to build the NLO versions of the BK [15] and B-JIMWLK [15][16][17][18][19] equations, and also of the "'impact factors" (the scattering matrix elements void of evolution) for a few, relatively simple, processes, like single inclusive hadron production in pA collisions [20][21][22], the structure functions for electron-proton (ep) or electron-nucleus (eA) deep inelastic scattering (DIS) at small Bjorken x [23][24][25], exclusive diffractive vector-meson and dijet production in DIS at small x [26][27][28], and inclusive photon + dijet production in eA collisions at small x [29].
The original NLO results for both the BK equation and (some of) the impact factors turned out to be problematic, in the sense of generating instabilities in the high energy evolution [30], or negative cross-sections for particle production [31][32][33][34]. These difficulties have been eventually understood and cured. In the case of the BK equation, the solution [35][36][37][38][39][40][41][42] involves all-order resummations which enforce the proper time-ordering for the lifetimes of the soft gluon fluctuations of the dilute projectile (a parton from the proton, or a colour dipole in the case of DIS). The scheme-dependence of such resummations is considerably reduced when using the rapidity of the dense target (the large nucleus) as the "evolution time" [40]. For the NLO impact factors, the negativity problem arises when enforcing a separation between the leading-order BK evolution and the NLO impact factor which is local in rapidity [43] (thus following the traditional prescription of the k T -factorisation [44,45]). So, the simplest way to avoid the problem -and obtain a positive-definite cross-section -is to give up this separation, that is, to keep the high-energy evolution and the NLO corrections together, in a generalised impact factor [43,[46][47][48]. In this "unsubtracted" scheme, the only evolution which needs to be factorised from the "hard" impact factor is the DGLAP evolution on the external lines -that is, the evolution of the parton distribution for the parton from the proton which initiates the scattering, and that of the fragmentation functions for the produced hadrons.
One process that has attracted much interest over the last years is the production of a pair of jets or hadrons in "dilute-dense" -proton(deuteron)-nucleus (pA) [49][50][51][52][53][54][55][56][57][58] or electron-nucleus (eA) [51,[59][60][61][62][63][64][65] -collisions at forward rapidities. This process could be used to probe saturation even in the (experimentally more accessible) set-up where the final jets/hadrons have transverse momenta significantly larger than the saturation momentum in the dense target. Indeed, the multiple scattering off the saturated gluons in the target is one of the mechanisms responsible for the transverse momentum imbalance between the final jets, hence for their distribution in the relative azimuthal angle. (The other important such a mechanism is the final state radiation, leading to the so-called "Sudakov factor" [66].) In turn, this imbalance leads to a broadening of the final particles distribution in the relative azimuthal angle ∆φ, around the back-to-back peak at ∆φ = π. Such a broadening has indeed been observed in d+Au collisions at RHIC [67,68], although at this level it looks difficult to distinguish the effects of saturation from those of the final state radiation [69].
That said, it would be important to have a more precise calculation of the cross-section for dijets (or dihadrons) production within the CGC formalism. This is our main purpose in this and a subsequent paper, where we will compute the respective NLO corrections to the impact factor for the case of pA collisions. Specifically, in this paper we shall present the "real" NLO corrections -those associated with a final state which involves three partons, out of which only two are measured. By integrating out the kinematics of the unmeasured parton, one generates a NLO correction to the cross-section for dijet production. In the companion paper, we will compute the corresponding "virtual" correctionsthose associated with one-loop corrections to the amplitude.
For simplicity and to avoid a proliferation of cases, we shall restrict ourselves to the quark channel, that is, we shall only consider processes in which the final state is obtained via radiation from an initial quark from the wavefunction of the incoming proton. The quark channel is expected to dominate the forward particle production in the kinematics at RHIC, but the gluon channel should be added too in view of realistic applications to the LHC. We plan to do so in a later paper.
The starting point for computing the "real" NLO corrections to the dijet cross-section is the leading-order (tree-level) cross-section for "trijet" (three partons) production, that we obtained in a previous paper [1], by using the formalism of the light-cone wavefunction (LCWF). So, for our present purposes, it should be enough to "integrate out" one (any) of the three final partons in our results in [1]. However, this is not what we shall do in practice. Indeed, the LO trijet cross-section is by itself a very complex quantity and our respective results in [1] were presented in a cumbersome way, which is not convenient for the present purposes. (We have realised that when trying to match our results for "real" and "virtual" corrections to dijet production, e.g. in order to check the cancellation of the infrared divergences.) So, our first step in this paper will be to re-derive the results for 3-parton production in a streamlined way and with a better strategy for organising the final expressions. In this process, we shall also correct several errors which occurred in our original results [1], that we were able to identify thanks to various tests to be discussed below.
In view of the above, the results for the LO trijet cross-section to be presented in this paper should be viewed as our final respective results, in replacement of the previous ones in [1]. These results are still formal, as they involve undone Fourier transforms in the transverse plane: the calculation of the LCWF is performed using the transverse coordinate representation, to take profit of the eikonal approximation. Hence, our expressions for the cross-section involve Fourier transforms relating the transverse coordinates of the three final partons (in both the amplitude and the complex conjugate amplitude) to their respective momenta, as measured in the final state. To actually compute these Fourier transforms -say, in view of applications to the phenomenology -one would also need explicit results for the partonic S-matrices probed by these particular processes. The techniques for computing such S-matrices -from numerical solutions to the BK and JIMWLK equations, or via mean field approximations to the latter -are well documented in the literature, but their discussion goes beyond our present purposes. In agreement with previous studies [51,70], we shall find that in the multicolour limit N c → ∞, the various S-matrices that we shall encounter are built with just two non-trivial colour structures: the dipole and the quadrupole.
Given these results for trijet production, the "real" NLO corrections to dijets are easily obtained (at least, formally) by integrating over the longitudinal and the transverse momentum of one of the three final partons and then convoluting with the proton distribution function for the incoming quark. If one measures two hadrons in the final state, convolutions with the appropriate parton-to-hadron fragmentation functions are also needed.
As a check of our calculations, we shall study two special limits for which the results are a priori known: the case where the unmeasured parton is a soft gluon and that where this parton is produced via a collinear splitting. In the first case, we expect to recover (one step in) the B-JIMWLK evolution 1 of the LO cross-section for producing a quark-gluon pair. In the second case, one should find the DGLAP evolution for the incoming quark distribution and for the fragmentation functions of the produced quark and gluon. We indeed recover these expected results, but only after rather non-trivial manipulations, which in particular involve cancellations between many terms. Thus the good results that we find in these limits provides a rather stringent test on our results.
We shall not insist in separating out the LO evolution from the NLO impact factor. In the case of the B-JIMWLK evolution, this is in line with our preference for the "unsubtracted" scheme [43] alluded to above, which avoids potential problems associated with a local subtraction in rapidity. In the case of the DGLAP evolution, the subtraction of the collinear divergences is indeed mandatory, but it requires more refined techniques like the dimensional regularisation, that we shall develop in the companion paper devoted to "virtual" NLO corrections.
This paper is structured as follows. In Sect. 2 we shall briefly recall the result for the LO dijet cross-section in the quark channel, i.e. the Born-level cross-section for producing a quark-gluon (qg) pair at forward rapidities. This result is well known [1,49,51] but we shall often need it for comparisons with the NLO results to be obtained later. In Sect. 3, we present the LO (tree-level) results for the production of three partons: qqq and qgg. We first show the three-parton components of the quark LCWF in the final state and then use them to compute the relevant cross-sections. As already mentioned, in doing that we shall both reorganise and correct our original results in [1]. In Sect. 4 we compute the "real" NLO corrections to dijet production in the quark-initiated channel by (formally) "integrating out" one of the three partons in the trijet results. We also explain the simplifications which occur in the colour structure (i.e. in the partonic S-matrices) due to the fact that one of the partons is not measured. In Sect. 5 we show that in the limit where the unmeasured parton is a soft gluon, our NLO results reproduce, as expected, the (real part of the) JIMWLK evolution for the qg dijet cross-section. Finally, in Sect. 6, we consider the limit where the unmeasured parton is produced by a collinear splitting. We show that, in this limit, our NLO results develop collinear singularities, that we isolate to leading logarithmic accuracy and verify that they can be interpreted as one-step in the DGLAP evolution of the initial (q) and final (qg) states.

Dijet cross-section at leading-order
As a warm-up, let us briefly recall some steps in the derivation of the leading-order (LO) result for the cross-section for dijet production in pA collisions at forward rapidity (see e.g. [1,49] for more details). As mentioned in the Introduction, we consider only the quark channel: the parton collinear with the incoming proton and which initiates the process is taken to be a quark. (The corresponding result for the gluon channel can be found in Refs. [1,53].)

Kinematics
We start with some generalities on the kinematics and use this opportunity to introduce some of the notations. We chose the z direction along the collision axis and work in a frame in which the proton is an energetic right-mover with (light-cone) longitudinal momentum Q + , whereas the nucleus is an ultrarelativistic left mover, with longitudinal momentum P − per nucleon. (We shall neglect the nucleon masses in what follows.) We more precisely assume that the nucleus target carries most of the total energy, so the high-energy evolution via the successive emissions of soft gluons is fully encoded in the nuclear gluon distribution.
To leading order, dijet production at forward rapidities and in the quark channel proceeds as follows (see also Fig. 1): a quark initially collinear with the proton, with longitudinal momentum q + = x p Q + , scatters off the dense gluon system in the nuclear target and emits a gluon in the process. The two "jets" are the final quark, with longitudinal momentum p + = x 1 Q + and transverse momentum p, and the emitted gluon, with k + = x 2 Q + and final transverse momentum k. The target has zero "plus" momentum, hence the respective component is preserved by the scattering: x p = x 1 + x 2 ; accordingly, we shall also write k + = ϑq + and hence p + = (1 − ϑ)q + , with ϑ = x 2 /x p ≤ 1 the gluon splitting fraction. On the other hand, the collision can transmit a transverse momentum of the order of the saturation momentum Q s -the typical momentum of a gluon at saturation -, hence we expect an imbalance |p + k| ∼ Q s between the final jets.
The two jets are put on-shell by the collision, hence they must receive from the nucleus a "minus" component equal to their total light-cone energy; writing this as a fraction x g of P − , we deduce with s = 2Q + P − (the center-of-mass energy squared of the collision). It is customary to express the longitudinal fractions x p and x g in terms of the (pseudo)rapidities η 1 ≡ (1/2) ln(p + /p − ) and η 2 ≡ (1/2) ln(k + /k − ) of the produced jets in the center-of-mass frame, where Q + = P − = s/2. Using p + = (p ⊥ / √ 2)e η 1 and similarly k + = (k ⊥ / √ 2)e η 2 (with p ⊥ ≡ |p| etc), one finds The forward dijet kinematics corresponds to the situation where η 1 and η 2 are both positive and larger than 1. In this regime, one has x g x p < 1, showing that the forward particle production explores the small-x g part of the nuclear wavefunction. The relevant value of the nuclear saturation momentum increases with decreasing x g , due to the rise in the gluon density in the transverse plane, via soft gluon emissions. This is governed by the non-linear evolution equations (BK or B-JIMWLK), which imply , with λ s 0.20 (see e.g. [41] for a recent study). Note that in deriving Eq. (2.2) we did not necessarily assume 2 → 2 kinematics (i.e. qg → qg). In fact, the most interesting case for us here is that of multiple scattering, which involves the exchange of arbitrarily many soft gluons between the quark-gluon fluctuation of the proton and the nucleus. In such a case, x g P − is the total LC energy transmitted from the target to the quark-gluon system and p + k is similarly the total transverse momentum.
Under the present assumptions, multiple scattering can be resummed to all orders within the eikonal approximation: the nuclear target appears to the proton as a Lorentz-contracted shockwave (say, localised at x + = 0) and the transverse coordinates of a projectile parton (quark or gluon) is not modified by the collision. It is then convenient to work in the transverse coordinate representation and Fourier transform to transverse momenta only at the end of the calculation. In this representation, the only effect of the scattering is a colour precession of the parton wavefunction, represented by a Wilson line in the appropriate representation of the colour group SU(N c ).

The quark-gluon component of the quark light-cone wavefunction
We use the light-cone wavefunction (LCWF) formalism and the projectile LC gauge A + = 0. The initial state at (LC) time x + → −∞ is taken to be a bare quark: |q in = |q α λ (q + , q) , with q + , q, λ and α denoting the quark longitudinal and transverse momenta, its spin, and its colour state, respectively. (Our conventions for the bare Fock states and for the action of the associated creation and annihilation operators are presented in Appendix A.) The outgoing state at x + → ∞ is computed as |q out = U I (∞, 0)Ŝ U I (0, −∞) |q in , where the QCD evolution operators U I (0, −∞) and U I (∞, 0) describe QCD radiation prior and after the scattering, respectively, andŜ is the S-matrix operator for the scattering between the parton system which exists at time x + = 0 and the shockwave.
To compute quark-gluon production to leading order (LO), it is enough to expand the action of the evolution operators to O(g), e.g.
where H int is the interaction part of the Hamiltonian, |i , |j are energy-momentum eigenstates of the free QCD Hamiltonian H 0 (i.e. the Fock states built with bare partons) and E i , E j are the corresponding LC energies (the sum of the minus components of the partonic 4-momenta for all the partons composing a Fock state, assumed to be on-shell). Taking |i = |q α λ (q + , q) , the outgoing state is a superposition of a single-quark state (whose colour has been rotated by the scattering) and the quark-gluon state, in which we are primarily interested. The qg Fock-space component is most convenient written in the transverse coordinate representation, where the matrix elements ofŜ are diagonal; it then reads [1]: where w is the transverse coordinate of the incoming quark (conjugated to its transverse momentum q), x and y are the transverse coordinates of the final quark and gluon, and the constraint w = (1−ϑ)x+ϑy follows from the conservation of the transverse momentum. Furthermore, g b i (ϑq + , y) is a gluon with transverse helicity i, colour state b, longitudinal momentum ϑq + and transverse position y. We also used the following notation for the spinor/helicity structure: The two terms inside the square bracket in Eq. (2.4) come from the action of the S-matrixŜ and refer to the situation when the scattering with the shockwave occurs before and respectively after the gluon emission (see Fig. 1). U ba (x) and V βα (x) are Wilson lines in the adjoint (for the gluon) and respectively fundamental (for the quark) representation. In matrix notations, with the colour field A − a representing the small-x gluons in the target. In the CGC effective theory, this field is random and must be averaged out at the level of the cross-section.
The outgoing state in Eq. (2.4) is clearly vanishing in the absence of scattering, i.e. in the limit where U → 1 and V → 1. Indeed, in the absence of any interaction (like the scattering off a nuclear target), the quark-gluon pair cannot be produced in the final state (since an on-shell quark cannot emit a gluon). Furthermore, using the identity one sees that the difference of Wilson lines in Eq. (2.4) also vanishes when y → x. This is the expression of colour conservation: in the limit where the transverse separation between the daughter partons shrinks to zero, the quark-gluon pair scatters off the nuclear target in the same way as its parent quark, hence the bare quark and its qg fluctuation cannot be disentangled by the collision.

The quark-gluon cross-section at leading order
The inclusive cross section for forward quark-gluon production in quark-nucleus (qA) scattering can be computed as the average of the respective product of number density operators over the momentumspace version of the outgoing state, as obtained from Eq. (2.4) via a Fourier transform: (2.8)  Specifically, with our present conventions one can write As already mentioned, we can choose q = 0 (the longitudinal axis is taken to be parallel to the collision axis). After some simple algebra, one finds the following result for the forward quark-gluon production cross section at LO [1,49]: with C F ≡ t a t a = N 2 c −1 2Nc and a compact notation for the transverse integrations: x ≡ d 2 x. It is understood that ϑ = k + /q + , w = (1 − ϑ)x + ϑy and w = (1 − ϑ)x + ϑy. To obtain this result, we have also used the identity ε il ε jk = δ ij δ kl − δ ik δ jl to perform the sums over the quark helicities and the gluon polarisations (recall Eq. (2.5)): In writing Eq. (2.10), we have also performed the average over the random colour fields A − describing the soft gluons from the target and we have introduced the following S-matrices describing the forward scattering of colourless systems made with up to four partons: a quark-antiquark dipole, a quark-antiquark-gluon triplet (this is illustrated in Fig. 2), and finally a quark-antiquark pair accompanied by two gluons: Q(x, y, y, x) S(y, y) . (2.14) The second equalities in the r.h.s. of Eqs. (2.14) and (2.13) are obtained after using Eq. (2.7) together with standard Fierz identities. Besides the colour dipole already introduced in Eq. (2.12) they also involve the quadrupole, The final, approximate, equalities in Eqs. (2.14) and (2.13) hold in the multi-colour limit N c → ∞, which allows for important simplifications, as already visible in the above results. It is understood that the target averages occurring in the above equations (2.12)-(2.15) are computed over the nuclear gluon distribution evolved down to a longitudinal momentum fraction x g , cf. Eq. (2.1). So, in principle, these S-matrices must be obtained from solutions to the non-linear B-JIMWLK equations, which in general form an infinite hierarchy. At large N c , one can use simpler equations, which are closed: the BK equation for the dipole S-matrix together with the evolution equation for the quadrupole obtained in [71][72][73]. The latter is still cumbersome to use in practice, so it is useful to notice that a mean field approximation relating the quadrupole to the dipole [51,54,[72][73][74] appears to work quite well, even for finite N c [75].
Given the partonic cross-section (2.10), the contribution of the quark channel to the leading-order cross-section for dijet production in pA collisions at forward rapidity is obtained as where q f (x p , µ 2 ) is the quark distribution function of the proton, for a longitudinal momentum fraction x p = q + /Q + and a resolution scale µ 2 . At LO, the value of x p is fixed by the δ-function in Eq. (2.10). Furthermore, µ 2 should be chosen of the order of the hardest among the transverse momenta, k 2 or p 2 , of the produced jets. The final result at LO and for large N c reads where q + = p + + k + , ϑ = k + /q + , and x p = q + /Q + . We shall later need also the expression of the LO cross-section for the case where one is measuring a pair of hadrons (instead of jets) in the final state. Under the assumption that the parton → hadron splitting is collinear, this cross-section is obtained by convoluting Eq. (2.17) with the quark → hadron and gluon → hadron fragmentation functions: where it is understood that p 1 = p/z 1 , k 1 = k/z 2 (for both longitudinal and transverse components), x p = (p + 1 + k + 1 )/Q + , and ϑ = k + 1 /(p + 1 + k + 1 ). This expression (2.18) for the dijet cross-section is a particular example of the hybrid factorisation, originally introduced in the context of single hadron production in pA collisions [2][3][4]. Thus is "hybrid" in the sense that the initial and final states are treated in the spirit of the collinear factorisation, whereas the hard process is rather computed using the CGC effective theory for QCD at high energy (which in turn can be viewed as a generalisation of the k T -factorisation originally developed in the context of the BFKL evolution [44,45]).

Trijet cross-section at leading order
In preparation for the calculation of the real NLO corrections to dijet production, let us first revisit our results for the LO trijet (qqq and qgg) cross-sections, originally presented in Ref. [1]. As compared to [1], we shall improve these results at two levels: (i) we shall rewrite them in a different way, which is not only more compact and physically more transparent, but also better suited for the purposes of the NLO calculation, and (ii) we shall correct several errors and misprints, that we failed to identify in Ref. [1], by lack of appropriate tests. Hence, the results for the trijet cross-section to be presented in this section should be viewed as our final respective results, in replacement of those in Ref. [1]. And as matter of facts, we have obtained these new results via an independent calculation, and not by just reshuffling terms in the original results from Ref. [1].

The tri-parton components of the quark outgoing state
We first present the expressions of the tri-parton Fock-space components of the quark outgoing state. These are obtained by expanding the QCD evolution operators in |q out = U I (∞, 0)Ŝ U I (0, −∞) |q in to second order in the interaction Hamiltonian (see e.g. Ref. [1] for details). In this expansion, we shall ignore the virtual corrections for the time being (they will be computed in a subsequent paper); that is, we shall only keep those second-order terms which contribute to a tri-parton (qqq and qgg) final state. We shall not present the details of the calculations -they would be very similar to those described at length in Ref. [1] -but only emphasise the differences w.r.t. [1], which mainly refer to a reshuffling of terms. For the convenience of the reader, our conventions for the light-cone wavefunction formalism are summarised in Appendix A, whereas Appendix B exhibits all the matrix elements of the QCD interaction Hamiltonian which are relevant for this calculation.

The tri-quark final state
To explain this reshuffling, we shall focus on the three-quark (qqq) final state. As explained in [1], this state receives two types of contributions, 2) The three terms in the r.h.s. of this equation correspond to the different possible insertions of the shockwave relative to the two parton branchings, as illustrated by the Feynman graphs in Fig. 3. In the first term, cf. Fig. 3.a, both parton branchings occur prior to the scattering off the shockwave, hence they are generated by the second-order expansion of the evolution operator U I (0, −∞). Accordingly, both energy denominators involve differences w.r.t. the energy E q of the initial state. Similarly, the second term, cf. Fig. 3.b, describes the process where both branchings occur in the final state (after the scattering), hence they are generated by U I (∞, 0). The corresponding energy denominators involve differences w.r.t. the energy E qqq of the final state. Finally, the third term, cf. Fig. 3.c, describes one emission prior to the scattering, as generated by the first order term in U I (0, −∞), and a second one after the scattering -the first order term in U I (∞, 0). The reshuffling of terms is based on the following identity: which has a simple physical meaning: it shows that in the absence of scattering (i.e. in the limit S → 1), the three second-order corrections in Eq. (3.2) exactly cancel each other. In turn, this is a consequence of the fact that an on-shell quark cannot radiate gluons in the absence of any scattering  (see the discussion in [1]). In what follows, we shall rely on this identity to re-express the middle term in Eq. (3.2), cf. Fig. 3.b, as the sum of two negative contributions which exhibit the same energy denominators as the two other terms there (see Fig. 4 for a graphical illustration of this manipulation). After this rewriting, the contribution of the initial-state scattering is treated as a subtraction term for the two other contributions (those shown in Figs. 3.a and 3.c).
It is in principle straightforward to implement this reshuffling of terms via manipulations of the original results for the trijet final state, as presented in [1]. However, in presenting our respective results in what follows, we shall often use different integration variables compared to [1]. Moreover, as stated at the beginning of this section, our original results in [1] contained some misprints and errors, that we have identified (and corrected) after explicitly working out special limits. In view of that, it is perhaps not so useful to compare with Ref. [1] anymore, but simply rely on the new results to be presented here. This is in particular true for the qqq component of the outgoing quark state. In Ref. [1], this has been obtained by separately evaluating the three terms in the r.h.s. of Eq. (3.2), with the result shown in Eq. (4.9) of that paper. After treating the second term there (which corresponds to Fig. 3.b) as a subtraction term, performing some changes of variables, and correcting an error, one finds Our notations are illustrated in Fig. 3): ϑ and ξ are the longitudinal momentum fractions (w.r.t. the momentum q + of the incoming quark) of the intermediate gluon and of the quark produced by the splitting g → qq, respectively. Also, x and y are the transverse coordinates of the quark and the emitted gluon, whereas z and z similarly refer to the quark and the antiquark generated by the gluon. Notice that y is not an independent coordinate, since it must coincide with the center of energy of the qq pair produced by the gluon decay. Similarly, the transverse coordinate w of the incoming quark must be the same as the center of energy of the quark-gluon pair in the intermediate state, or of the three quarks in the final state. Explicitly, The δ-function inside the integrand enforces the condition on w shown in the second equality in Eq. (3.5), thus introducing a constraint on the 5 integrations visible in the r.h.s. of Eq. (3.4). The transverse coordinates denoted with capital letters, that is, are the transverse separations between the daughter partons after each of the emission vertices. Finally, the spinorial structure of the two branching vertices is encoded in the following tensors: The two terms within the square brackets in Eq. (3.4) represent the contributions of the two graphs in Fig. 3.a and Fig. 3.c respectively, minus the respective subtraction terms generated by the graph in Fig. 3.b, as shown in Fig. 4. Thanks to these subtractions, the numerator in each of these terms vanishes in the limit where the transverse separations between the daughter partons (as appearing in the respective denominator) are shrinking to zero. Consider e.g. the second term: when R = x−y → 0, meaning that the three coordinates x, y and w become coincident with each other, one can use the following identity (cf. Eq. (2.7)) to verify that the respective numerator vanishes, as announced. Similarly, the numerator in the first line vanishes when both R → 0 and Z → 0 (meaning that all transverse coordinates become degenerate: w = x = y = z = z ). Furthermore the two terms within the square brackets in Eq. (3.4) cancel each other in the limit Z → 0 (i.e. z = z = y) at fixed R. To see this, one should use yet another version of the identity Eq. (2.7), namely, which expresses the fact that the scattering cannot distinguish a qq fluctuation of zero size from its parent gluon (recall the discussion after Eq. (2.7)). These properties of the LCWF will be later useful in demonstrating the cancellation of "ultraviolet" (short-distance) singularities in the calculation of cross-sections for particle production. Such properties are more difficult to check on the original expression for the LCWF, i.e. Eq. (4.9) in Ref. [1].
The previous discussion also implies that the second term in Eq. (3.4) can be generated from the first term there via the simultaneous replacements z → y and z → y. This allows us to introduce a simpler notation, that we shall systematically use throughout this paper. Specifically, Eq. (3.4) can be equivalently written as (3.10) Besides this "regular" contribution, which involves an intermediate state with a propagating gluon, there is an additional contribution built with the instantaneous piece of the gluon propagator (see Fig. 5). This is given by Eq. (4.13) in Ref. [1], which with the present conventions reads where c ≡ (1 − ϑ)x + ξz + (ϑ − ξ)z is a compact notation for the center of energy of the final quark triplet. Note that, in this case, there is no intermediate gluon state: the cut gluon propagator in Fig. 5 together with the two attached QCD vertices defines a local effective vertex for the decay q → qqq of the initial quark into the three final ones. So, in particular, there is no analog of the second term in Eq. (3.4) (the one denoted as (z, z → y) in Eq. (3.10)). That said, it is still convenient to introduce a compact notation (similar to that in Eq. (3.10)) for the complete triquark outgoing state, i.e. the sum of the two contributions in Eqs. (3.4) and (3.11). This reads is an effective vertex for the splitting of the original quark into the three final quarks, which includes contributions from both propagating and instantaneous intermediate gluons: It is understood that the second piece in Eq. (3.13) vanishes after substituting z → y and z → y : Hence the subtraction of the contribution from the intermediate gluon state is truly performed only for the case of a propagating gluon, as it should. Ultimately, one should keep in mind that the compact expression for the three-quark outgoing state in Eq. (3.12) is merely a convenient notation, which is rather formal, but will be useful to simplify the writing for our final results.

The final state with one quark and two gluons
In this subsection we shall exhibit the remaining 3-parton Fock space components of the quark outgoing LCWF: those in which the original quark is accompanied by two gluons. There are two possible topologies for the gluon emissions, as shown in Figs. 6 and 7, respectively. Besides these "regular" graphs, where the intermediate parton (quark or gluon) is propagating, there are similar graphs which involve the instantaneous piece of the respective propagator. (These are not shown, but can be easily inferred by comparing with Fig. 5.) All these graphs have been computed in [1], with results that we shall here rewrite by regrouping terms and also correct whenever necessary, as explained in the previous subsection on the example of the three-quark final state.
We start with the graphs where both gluons are emitted by the original quark, cf. Fig. 6. Their contribution is shown in Eq. (4.14) of [1], that can be rewritten as  17) are the transverse separations between the daughter partons after each emission vertex. The spinorial structure is encoded in φ il λ 1 λ (ξ) (cf. Eq. (3.7)) for the first emission vertex and respectively for the second vertex; notice the relation τ ij λ 1 λ (ϑ, 0) = φ ij λ 1 λ (ϑ). The structure of Eq. (3.15) is reminiscent of that in Eq. (3.10). The first term within the squared brackets represents the contribution of the final-state interaction, cf. Fig. 6.a, minus a piece generated by the initial-state interaction, cf. Fig. 6.b. The second term, as obtained from the first term via the replacements x, z → y, refers to the intermediate-state interaction, cf. Fig. 6.c, minus the remaining contribution of Fig. 6.b.
The contribution of the corresponding instantaneous graph, for which there is no analog of Fig. 6.c, is found as As in the previous subsection, it is convenient to group together "regular" and "instantaneous" contributions in a unique expression involved a non-local effective vertex -here, for the splitting of the original quark into a qgg system. Then the sum of Eq. (3.15) and Eq. (3.19) reads It is understood that the second piece in Eq. (3.21) vanishes when X = x − z → 0. The remaining topology for a final state with two gluons and a quark is very similar to that previously discussed for the tri-quark final state (compare the graphs in Figs. 7 and 3, respectively). Hence, we shall simply present here the final result -the analog of Eq. (3.12) -without further explanations. This reads where the first term in the r.h.s. is the three-gluon vertex whereas the second term comes from the graphs with an instantaneous intermediate gluon.
To obtain an explicit expression for the subtracted term denoted as (z, z → y) in Eq. (3.22), the following identity is useful: As before, one can check that, thanks to this subtraction, the two terms within the square brackets cancel each other in the limit Z → 0.

The trijet cross-section
The leading-order cross-section for the inclusive three partons ("trijet") production in the quark-nucleus collision is computed similarly to Eq. (2.9), that is, as the expectation value of the product of three number-density operators (themselves built with Fock space operators for bare partons).

The qqq final state
For the qqq final state, one writes (a factor 2πδ(k where as shown in the second line, only the qqq Fock-state component of the outgoing quark state, cf. Eq. (3.12), is involved in this calculation. A priori there are four possible contractions for the product N q (k 1 )N q (k 2 ) of quark number density operators. However, two of these contractions correspond to non-planar graphs -these are graphs where the original parton on one side of the cut is contracted with the quark generated by the gluon decay on the other side of the cut (see Fig. 13.b in [1]) -which are suppressed in the multicolour limit N c → ∞. In what follows we shall systematically work in this limit, since it allows to simplify the colour structure of our results and also to render them physically transparent. Hence we shall keep only two of the four possible contractions -those corresponding to planar graphs as illustrated in Fig. 8. We shall explicitly write the result for the case where the leading quark has momentum k 1 and the other quark has momentum k 2 . The other term can be simply obtained by permuting these two momenta.
A straightforward calculation using the outgoing state in Eq. (3.12) together with the Fock space rules summarised in Appendix A yields 2 (3.27) The notations here are similar to those in Eq. (3.12). The transverse coordinates x, z, z , y, w refer to the direct amplitude (DA) and have the same meaning as in Fig. 3 and Eq. (3.5). We use a bar to indicate the corresponding coordinates in the complex conjugate amplitude (CCA). We use the notations in Eq. (3.6) for the transverse separations between the emitted partons, that is, The longitudinal momentum fractions ϑ and ξ are the same in the DA and in the CCA, since they are fully fixed by the kinematics of the final state, as follows: The δ-function enforcing longitudinal momentum conservation implies k + 3 = (ϑ − ξ)q + , which in particular requires ϑ ≥ ξ (i.e. k + 1 + k + 2 ≤ q + ). Let us now explain the new structures occurring in Eq. (3.27). The tensorial kernel is defined as with the effective vertex Φ ij λ 3 λ 2 λ 1 λ introduced in Eq. (3.13) (the star in Φ * denotes complex conjugation). The product of effective vertices in the numerator can be explicitly computed as (3.31) As emphasised by our notation, in evaluating this product we have excluded the instantaneous pieces from the two effective vertices. These pieces have a very simple tensorial structure, so they can be easily inserted when needed 3 .
2 In the large-Nc limit under consideration, the quark Casimir factor apparent in Eq. (3.27) can be as well approximated as CF Nc/2; this applies to all the results involving CF to be shown in this paper. 3 These instantaneous pieces do not contribute to either the soft, or the collinear, limit that we shall study later; so, a result like Eq. (3.31) is indeed sufficient for our purposes in this paper. The effects of the collision are encoded in the function W 0 , defined as the following linear combination of partonic S-matrices: (3.32) Our notations for the partonic S-matrices are intended to describe (via the lower scripts) the partonic composition of the multi-parton system which scatters off the shockwave and to also distinguish (via a bar on the transverse coordinates) between partons in the DA and in the CCA, respectively. The ordering of the transverse coordinates in the argument follows that of the lower subscripts. The indices/arguments corresponding to partons in the DA appear on the left to those representing the CCA. When reading the lower indices, one should also keep in mind that a quark in the CCA formally counts like an antiquark (from the viewpoint of the direction it its colour charge flow). According to these conventions, the S-matrix S qqqqqq (x, z, z , x, z, z ) refers to the eikonal scattering of a system made with 6 quarks: 2 quarks (x, z) and one anti-quark (z ) in the DA, together with 2 "anti-quarks" (x, z) and one "quark" (z ) in the CCA. This describes the situation where the collision occurs in the final state (i.e. after the second splitting) in both the DA and in the CCA (see Fig. 8.left), and reads where the approximate equality holds for large N c : in this limit, the 6-quark S-matrix factorises into the product of a dipole, S(z, z), and a quadrupole, Q(x, z , z , x). Furthermore, S qqqq (x, z, z , w) and S qqqq (w, x, z, z ) represent interference terms where the collision with the shockwave occurs in the final state on one side of the cut, and in the initial state (prior to the first splitting) on the other side; e.g., (3.34) Finally, the colour dipole S (w, w) ≡ S qq (w, w) describes initial-state interactions in both the DA and the CCA. The normalisation factors in the above definitions are such that all the individual S-matrices reduce to unity in the absence of scattering. Accordingly, both the linear combination in Eq. (3.32) and the cross-section (3.27) vanish in that limit, as expected. The overall colour structure becomes remarkably simple at large N c : The three other terms within the square brackets in Eq. (3.27), which are obtained as various limits of the first term, refer to situations where the scattering in the final state gets replaced by scattering in the intermediate state, as explained in relation with Eq. (3.10). E.g. the second term within the square brackets, denoted as z, z → y , is illustrated in Fig. 8.right. We recall that this term can be obtained by letting Z → 0 in the kernel (3.30) (which in particular means using the simplified version of the effective vertex shown in Eq. (3.14)) and replacing z → y and z → y in the arguments of the partonic S-matrices in the r.h.s. of Eq. (3.32).
It is quite instructive to display the version of Eq. This S-matrix structure is identical to that occurring in the integrand of Eq. (2.17) for LO quark-gluon production (qA → qg +X). This is easy to understand: in both cases, the interactions with the nuclear shockwave refer to the partons involved in the branching q → qg. So far, we have not specified the x g -argument of the various S-matrices in Eq. (3.32), i.e. the "minus" longitudinal momentum fraction of the gluons from the nuclear target which are involved in the production of the three-parton state. Clearly, this is given by the generalisation of Eq. (2.1) to a final state involving three partons, that is, where x i = k + i /Q + are the "plus" longitudinal momentum fractions of the produced partons. Needless to say, this value for x g applies not only to the qqq final state discussed so far, but also to the qgg final states to which we now turn.

The qgg final state
For the qgg final state, the trijet cross-section is computed as (once again, a factor 2πδ(k + 1 +k + 2 +k + 3 −q + ) is implicitly understood in the l.h.s.) The second line involves the qgg Fock component of the quark outgoing state, which, as explained in Sect. 3.1.2, is the sum of two contributions, corresponding to two different topologies for the gluon emission vertices, as illustrated in Fig. 6 and Fig. 7, respectively. It is therefore natural to split the cross-section in Eq. (3.38) into three contributions. In the first one, the topology in Fig. 6 is used for the quark LCWF in both the DA and the CCA. Similarly, the second contributions involves the topology in Fig. 7 alone. Finally, the third contribution represents interferences between the two topologies. Besides, each of these three contributions is built with two pieces, corresponding to the two possible permutations for the momentum labels of the produced gluons. As before, we shall work in the limit of a large number of colours, in which one discards the non-planar graphs. One should however pay attention to the fact that the symmetry of the triple gluon vertex makes that the graphs generated by exchanging the two final gluons in Fig. 7 are still planar, and therefore must be kept in the final result. This will result in appropriate symmetry factors. A straightforward calculation using Eq. (3.20) yields the following expression for the first piece of the cross-section: is defined as with the effective vertex Ξ ijmn λ 1 λ introduced in Eq. (3.21) and the longitudinal momentum fractions ϑ and ξ shown in Eq. (3.29). The product of effective vertices in the numerator can be explicitly computed: Furthermore, W 1 denotes the following linear combination of partonic S-matrices: qqgg w, x, z, z + S (w, w) . (3.42) We use the same notations as explained after Eq. (3.32) in order to synthetically summarize the partonic content for both the DA and the CCA. The first S-matrix in the r.h.s., that is, describes final-state interactions in both the DA and the CCA and hence it includes Wilson lines for 6 partons: one quark (x) and two gluons (z, z ) in the DA, and one "anti-quark" (x) and two gluons (z, z ) in the CCA (see Fig. 9). Similarly, is an interference term between final-state and initial-state interactions. Once again, the large N c version of W 1 can be fully expressed in terms of colour dipoles and quadrupoles: As a simple check, we notice that after the double replacement x, z → y & x, z → y (which yields the last term within the square brackets in Eq. (3.39)), the r.h.s of Eq. (3.42) takes the same form as for the scattering of a quark-gluon pair 4 , as expected.
Consider also the second piece in the cross-section (3.39), as obtained by permuting the momenta k 2 and k 3 of the final gluons. This can be computed from the first piece (the one explicitly shown in Eq. (3.39) and illustrated in Fig. 9) by exchanging k 2 ↔ k 3 within the Fourier phases and k + 2 ↔ k + Figure 10. A graph contributing to the cross-section in (3.46); the associated S-matrix is shown in Eq. (3.50).
We now turn to the second contribution to the cross-section for qgg production, that involving the topology in Fig. 7. We now use Eq. (3.22) to deduce dσ qA→qgg+X whose numerator can be explicitly computed as (as before, we exclude the contribution of the instantaneous pieces of the vertices, to simplify writing) Furthermore, W 2 is the following linear combination of partonic S-matrices: qqgg w, x, z, z + S (w, w) , (3.49) Figure 11. A particular graph contributing to the "interference" cross-section (3.52); its colour structure is exhibited in Eq. (3.55).
where the 6-parton S-matrix describes final-state interactions in both the DA and the CCA, (the two terms correspond to permutations of the final gluons), whereas the 4-parton S-matrix describes the interference between scattering in the final state and the initial state, respectively: As a check, one can see that after performing the double replacement z, z → y & z, z → y in Eq. (3.49), one recovers the colour structure describing the scattering of a quark-gluon pair (the structure shown in Eq. (3.36) at large N c ). Note that, as compared to Eq. (3.39), Eq. (3.46) contains an additional factor of 4 (besides the modified colour factor, which reflects the different structure of the partonic S-matrices). This factor of 4 is related to the symmetry of the amplitude in Fig. 7 under the exchange of the two final gluons. This in turn implies that the integrand of Eq. (3.46) is symmetric under the simultaneous exchanges z ↔ z , z ↔ z , and ξ ↔ ϑ − ξ, as it can be easily verified. We used this symmetry property twice: (i) for a given assignment of the momenta of the two gluons (e.g., k + 2 = ξq + and k + 3 = (ϑ − ξ)q + , as in Fig. 10), there are two possible Wick contractions in the calculation of the expectation valuê N g (k 2 )N g (k 3 ), cf. Eq. (3.38), which give identical results; (ii) the graphs obtained by permuting the gluon momenta, k 2 ↔ k 3 , give identical results as well, because we can undo the effect of this permutation via a change of variables which leaves the integrand unchanged.
The remaining terms refer to interferences between the two topologies shown in Fig. 6 and Fig. 7, respectively. We shall explicitly show the piece where the topology in Fig. 6 counts for the DA and that in Fig. 7, for the CCA (see Fig. 11 for an illustration). The total result can be then obtained by taking the double of the real part of the shown contribution. One finds dσ qA→qgg+X (3.52) It is important to notice that, in the subtracted terms, y is the function of x and z defined in Eq. (3.16), whereas y is rather a function of z and z defined by the "bared" version of Eq. (3.5). Furthermore, The combination of Wilson lines describing the scattering has the same general structure, that is, with however some new ingredients, namely, It might be interesting to notice that the last term in Eq. (3.52), as obtained via the double replacement x, z → y & z, z → y , does not has the same colour structure as expected for a final quark-gluon state (unlike the respective limits for all the other colour functions W 0 , W 1 , and W 2 ). One finds indeed W 3 (y, z, y, x, y, y) S(y, z)S(y, x) S(z, y) − S(y, z) S(z, w) − S(w, y) S(y, x) + S (w, w) , which should be compared to Eq. (3.36). This difference is due to the fact that the pattern of the colour flow in the interference graphs is different as compared to the direct graphs. Note finally that, in the large-N c limit, the three functions encoding the colour structure in the various contributions to the cross-section, namely W 1 , W 2 , and W 3 , become very similar to each other: W 1 and W 3 take exactly the same form, shown in Eq. (3.45), whereas W 2 differs only through the additional symmetrisation with respect to the exchange of the two final gluons.
4 Next-to-leading order corrections: the real terms As mentioned in the Introduction, the next-to-leading order (NLO) corrections to the cross-section for forward dijet production in pA collisions can be divided into two classes: real and virtual. In the remaining part of this paper, we shall compute the real corrections -those associated with a final state which involves three partons ("jets"), out of which only two are measured. By integrating out the kinematics of the unmeasured parton, one generates a loop correction to the cross-section for the two measured jets. This loop opens in the direct amplitude (DA) and closes back in the complex conjugate amplitude (CCA). The virtual NLO corrections, on the other hand, refer to loop corrections to the amplitude itself. They will be addressed in a subsequent paper.

The di-quark jet production
We first consider final states which include two measured fermions: two quarks, or a quark-antiquark pair. The cross-section for di-quark jet production 5 is obtained by "integrating out" the final antiquark in our general formula for the three-quark final state, cf. Eq. (3.27): (4.1) (The subscript "rNLO" stays for real next-to-leading order corrections.) With reference to Eq. (3.27), it is quite clear that the integral over k + 3 can be trivially performed by using the δ-function for longitudinal momentum conservation, whereas the integral over k 3 yields a factor (2π) 2 δ (2) (z −z ), which allows one to identify the coordinates z and z of the unmeasured antiquark in the DA and the CCA, respectively. The result of Eq. (4.1) can be succinctly written as where it is understood that the trijet cross-section in the r.h.s. is given by Eq. (3.27), but without the δ-function expressing the conservation of longitudinal momentum. Some care must be taken, concerning the order of limits: the identification z = z must be made only after performing the subtractions which occur in the integrand of Eq. (3.27). For instance, in the subtracted term denoted as z, z → y , one must first perform the replacements z → y and z → y at fixed z and only then identify z with z . (The two limits do not commute with each other, as one can easily check.) The fact that the antiquark is not measured brings some simplifications in the structure of the 6-quark S-matrix in Eq. (3.33) (which describes final-state interactions, cf. Fig. 8-left): when z = z , the Wilson lines describing the scattering of the antiquark compensate each other by unitarity, V † (z )V (z ) = 1, and then the quadrupole appearing in the second line of Eq. (3.33) reduces to a dipole (we consider large N c , for simplicity): However, this simplification refers only to the first term within the square brackets of Eq. (3.27). Consider e.g. the second term, denoted as z, z → y . If one first replaces both z and z by y, and only then one identifies z = z , then the quadrupole structure survives in Eq. (3.33): S qqqqqq x, y, y, x, z, z = z O(x, y, z , x) S(y, z).  .29), which for the present purposes should be rewritten in terms of the longitudinal momentum fractions x 1 ≡ k + 1 /Q + and x 2 ≡ k + 2 /Q + of the measured partons ("jets") and the respective fraction x p ≡ q + /Q + of the original quark: As discussed in relation with Eq. (2.16), the physical dijet cross-section (to the order of interest) is obtained by averaging over x p with the quark distribution inside the proton: It is furthermore important to keep in mind that the cross-section (4.6) depends upon the "plus" longitudinal fractions x 1 , x 2 , and x p not only via its dependence upon ϑ and ξ, as explicit in Eq. (3.30) for the kernel K 0 , but also via the dependence of the various S-matrices upon the "rapidity" variable 6 x g , as introduced by their high-energy evolution. We recall that x g is the "minus" longitudinal momentum fraction carried by the gluons from the nucleus which participate in the scattering. This is indeed a function of x 1 , x 2 and x 3 = x p − x 1 − x 2 , as shown in Eq. (3.37).

The quark-antiquark di-jets
The case where the measured dijet is made with a quark-antiquark pair 7 is similarly obtained by integrating Eq. (3.27) over the kinematics of one of the two final quarks, e.g. over k 1 . This gives (4.7) 6 More precisely, the evolution "rapidity" related to xg is Yg ≡ ln(1/xg). 7 This case can be viewed as a NLO correction to the O(αs)-process in which a gluon collinear with the incoming parton splits into a qq pair while scattering off the nuclear target [1,53]. The respective cross-section involves the gluon distribution, which is however suppressed at large x, i.e. for the forward kinematics of interest for us here.
The above r.h.s. involves two contributions: the unmeasured quark k 1 can be either the incoming quark, or the quark produced by the decay of the intermediate gluon. The contribution of the first case to the qqq cross-section is explicitly shown in Eq. (3.27), while that of the second case can be obtained by exchanging k 1 and k 2 within the first contribution. Specifically, in the first case (k 1 corresponds to the leading quark), the integral over k 1 enforces x = x and the longitudinal momentum fractions should be evaluated as The fact that x = x simplifies the colour structure of the cross-sections in all the situations where the leading quark interacts in both the DA and the CCA: the respective Wilson lines compensate each other by unitarity. In what follows we shall exhibit these simplifications directly in the large N c limit. Specifically, the 6-quark S-matrix S qqqqqq (3.33) reduces to a product of two dipoles (one for each of the two measured partons): A similar simplification occurs in all the four terms within the squared brackets in Eq. (3.27); e.g., for the last term, as obtained via the double replacement z, z → y & z, z → y , one finds S qqqqqq (x, y, y, x = x, y, y) S(y, y) S(y, y). This is recognised as the large-N c version of the S-matrix for a colour dipole made with two gluons.
In the second case, where k 1 refers to the quark originating from the gluon decay, the integral over k 1 enforces z = z and the longitudinal momentum fractions are evaluated as As in the case of Eq. (4.2), the limit z → z does not commute with the various replacements in Eq. (3.27), like z, z → y , and must be performed after them. Accordingly, there is only one simplification in the colour structure, corresponding to the unique topology (in terms of shockwave insertions) in which the unmeasured gluon k 1 scatters both in the DA and in the CCA. One then has S qqqqqq x, z, z , x, z = z, z Q(x, z , z , x) . (4.12)

The two-gluon di-jets
The cross-section for the final state involving two measured gluons is obtained by "integrating out" the final quark in the trijet cross-section corresponding to the qgg final state. As discussed in Sect. 3.2.2, there are three contributions to the qgg cross-section, illustrated in Figs. 9, 10, and 11. By integrating out these three contributions over the kinematics (k + 1 , k 1 ) of the unmeasured quark, one finds (4.13) We again used the compact but somewhat formal notation introduced in Eq. (4.2), where it is understood that the δ-function expressing longitudinal momentum conservation is excluded from the trijet cross-sections. In turn, each of the three terms in the r.h.s. is made with two pieces, corresponding to the permutations of the measured gluons k 2 and k 3 . It is further instructive to notice the simplifications in the colour structure of the trijet crosssections which appear due to the fact that the quark is not measured (x = x). These simplifications refer to the S-matrices denoted as S (i) qggqgg with i = 1, 2, 3, which describe situations where the final quark scatters both in the DA and in the CCA. As usual, we show the ensuing simplifications only at large N c . In the three cases, one of the two quadrupole factors reduces to a dipole. One finds 8 , (4.14) and respectively

The quark-gluon dijets
The quark-gluon dijet final state already exists at leading order, as discussed in Sect. 2. The corresponding real NLO corrections are obtained by integrating out one of the two final gluons from the qgg trijet cross-section computed in Sect. 3.2.2, say the one with momentum k 3 . The complete result is therefore the sum of three terms, In turn, each of these three terms is made with two pieces, since the unmeasured gluon with momentum k 3 can be any of the two final gluons in Figs. 9, 10, and 11. Consider for definiteness the case where this is the gluon with longitudinal momentum ϑ − ξ and transverse coordinates z and z in the DA and the CCA, respectively. Then the integral over k 3 enforces z = z and the longitudinal fractions should be computed as in Eq. (4.11), that is, It is again instructive to display the main simplifications in the colour structure due to the fact that one of the final gluons is not measured (z = z ). One finds

(4.19)
8 As implied by its notation, the colour dipole S (z , z ) is built with a quark located at z in the CCA and an anti-quark located at z in the DA; recall Eq. (2.12).

Soft gluon emissions: recovering the B-JIMWLK evolution
In this section we shall study the behaviour of the previous results in the limit where one of the two final gluons is very soft, that is, either ξ → 0 or ϑ − ξ → 0. This is interesting since, in this limit, our general results should reduce to one-step in the B-JIMWLK (or BK) evolution of the leading-order dijet cross-section. As we shall see, this expectation is indeed verified and it provides a rather strong test of the correctness of our previous calculation. For more generality, we shall first address the soft gluon limit for the trijet (qgg) cross-section, before eventually specialising to the (real) NLO corrections to the dijet (qg) cross-section.

Direct emissions by the quark
We start with the contribution to the cross-section shown in Eq. (3.39), where the two gluons are directly emitted by the quark, both in the DA and in the CCA (recall Fig. 9). When (at least) one of these gluons is soft, there are two types of simplifications. First one can neglect the recoil of the emitter (here, the quark) at the respective emission vertex, meaning that its transverse coordinate is not modified by the soft emission. With reference to Fig. 6, when ξ → 0 one can approximate y w (1 − ϑ)x + ϑz , whereas for ϑ − ξ → 0, one has y x and w (1 − ϑ)x + ϑz.
Second, one can simplify the dependence of the kernel (3.40) -which, we recall, encodes information about both the emission vertices and the energy denominators -upon the longitudinal momentum fractions ϑ and ξ. It is easy to check that the instantaneous piece of the effective vertex Ξ (cf. Eq. (3.21)) does not contribute to the would-be singularity in the soft limit. Hence, one can use the expression (3.41) for the square of the effective vertex. This gives the following soft limits: Concerning the denominator of K 1 , the limit ϑ − ξ → 0 poses no special difficulty and yields The other eikonal limit, namely ξ → 0, is a bit more subtle, as it does not commute with the equal transverse points limits which define three of the four terms within the square brackets in Eq. (3.39).
Keeping the dominant terms as ξ → 0 in each of the respective kernels, one finds , when X = 0, but X = 0 , , when X = X = 0 , (5.3) These expressions exhibit single poles at ϑ = ξ and respectively ξ = 0, which are the expected infrared singularities associated with the bremsstrahlung of very soft gluons. What is however quite remarkable and also important for what follows, is the special way how these singularities appear within the structure of the cross-section (3.39): (i) the pole at ϑ = ξ is common to all the four terms in Eq. (3.39) and, moreover, the respective kernels become degenerate in this limit; (ii) the pole at ξ = 0 refers only to the last term in Eq. (3.39), as obtained via the double replacement x, z → y & x, z → y . These properties are essential in order to recover the JIMWLK evolution, as we now explain.
The limit ξ → 0 corresponds to the case where the soft gluon is the first one to be emitted (see Fig. 6). Property (ii) above tells us that the dominant graphs in this limit are the four graphs shown in Fig. 12, where the interaction with the shockwave occurs either in the initial, or in the intermediate, state (the soft gluon is drawn in red). Notice that this is only a small subset of all the possible topologies which contribute to the cross-section (3.39) in general: 4 out of 16. (Two of the 12 excluded diagrams are illustrated in Fig. 13.) What is special about the four graphs in Fig. 12 is the fact that there is no intermediate emission vertex between the shockwave and the emission of the soft gluon, neither in the DA, nor in the CCA (compare in this respect Fig. 12 to Fig. 13). In other terms, the parent quark radiates the soft gluon right before, or right after, its collision with the nucleus. This is in agreement with the fact that a soft gluon has a short formation time, hence it is first one to be emitted after a collision, or the last one to be emitted prior to it.
The graphs shown in Fig. 12 have the right structure to describe one step in the JIMWLK evolution of the LO dijet cross-section: they precisely encode the BK evolution of the last piece in Eq. (2.10) -the dipole S-matrix S (w, w) -, as we now demonstrate. Let us use the momentum assignments Figure 13. Diagrams which contribute to the qgg cross-section (3.39) in general, but which do not survive in the eikonal (soft) limit for the gluon with longitudinal momentum fraction ξ (shown in red). These graphs do not contribute to the JIMWLK evolution of any of the S-matrices which enter the LO cross-section (2.10).
shown in Fig. 9, i.e. k + 2 = ξq + and k + 3 = (ϑ − ξ)q + , with ξ → 0. In this limit, Eq. (3.39) simplifies to dσ qA→qgg+X where we have used w y (1 − ϑ)x + ϑz and we recall that Y = y − z w − z and X = x − z (and similarly for the respective coordinates with a bar). Also, we have ignored the soft momentum k + 2 in the δ-function for longitudinal momentum conservation. This is in the spirit of the soft gluon approximation to the high-energy evolution, which effectively violates energy conservation, albeit only marginally (by ignoring the soft gluons in the energy balance).
As expected, Eq. (5.4) is the same as the result of acting with the "production" version of the JIMWLK Hamiltonian [76][77][78] on the last term (the dipole S (w, w)) in the LO qg cross-section (2.10). The "production" Hamiltonian is a generalised version of the JIMWLK Hamiltonian which, when acting on the cross-section for particle production in dilute-dense (pA or eA) collisions, leads to the emission of an additional, soft, gluon, which is measured in the final state. If on the other hand this soft gluon is not measured, i.e. if one integrates out its kinematics, then one generates the "real" piece of the standard JIMWLK evolution. By integrating Eq. (5.4) over k + 2 = ξq + and k 2 , one finds Notice the simplifications in the colour structure due to the fact that the integral over k 2 has identified the coordinates z and z of the unmeasured gluon: z = z. The expression in the third line of Eq. (5.5) is recognised as the "real" part of the BK equation for the evolution of the dipole S (w, w). More precisely, it exhibits the integral version of the BK equation, that is, its formal solution as obtained by integrating the r.h.s. of that equation over ξ. This integral exhibits a logarithmic divergence at its lower limit (ξ = 0). In reality, this divergence is cut off by the condition of energy conservation. Indeed, this integral can be equivalently rewritten as (dξ/ξ) = (dx 2 /x 2 ), with x 2 = k + 2 /Q + ; then Eq. (3.37) together with the constraint x g ≤ 1 clearly implies a lower-limit on x 2 , hence on ξ.
This discussion also shows that the soft (ξ 1) part of the integral over the unmeasured gluon contributes to the high-energy (B-JIMWLK) evolution of the LO dijet cross-section, so it is only the remaining contribution at large ξ ∼ 1 which should be viewed as a genuine NLO correction to the "hard impact factor", i.e. to the partonic cross-section for the qA → qg + X process. That said, it turns out that the separation of the full integral over ξ between a "soft part" describing LO evolution and a "hard part" describing NLO corrections to the impact factor is generally subtle, due to the fact that the various S-matrices in Eq. (5.5) do implicitly depend upon ξ. (Indeed, their rapidity evolution scale is x g , which is a function of x 2 , hence of ξ, as shown in Eq. (3.37).) This dependence is essential, as argued in Refs. [43,47]: attempts to separate the "soft" evolution from the "hard" NLO impact factor which ignore this evolution may lead into troubles, like negative values for the NLO cross-section, as observed in the context of single inclusive particle production in both pA [21,31,46] and eA collisions [48]. A detailed discussion of such issues in the context of dijet production goes well beyond the purposes of this paper. But one should keep them in mind when trying to actually estimate the NLO corrections (say, for the purposes of the phenomenology).
Consider similarly the eikonal limit for the second emitted gluon, that with longitudinal momentum k + 3 = (ϑ−ξ)q + . Eq. (5.2) implies that, when this gluon is soft (ϑ−ξ 1), all the four terms within the squared brackets in Eq. (3.39) are multiplied by the same kernel, hence the respective colour structures can be simply added to each other. Remarkably, 12 of the 16 partonic S-matrices which are a priori Figure 15. Diagrams which contribute to the qgg cross-section (3.39) in general, but which do not survive in the eikonal (soft) limit for the gluon with longitudinal momentum fraction ϑ − ξ (shown in red).
contained in these structures mutually cancel in their sum, and we are left with dσ qA→qgg+X where X = x − z and X = x − z (we have used y x). The four surviving S-matrices correspond to the situations where the shockwave is inserted either in the final qgg state, or in the intermediate qg state (see Fig. 14), so that at least one gluon participates in the scattering (together with the quark, of course). For more clarity, we also show in Fig. 15 two graphs whose contributions have cancelled in the overall sum. The four surviving graphs in Fig. 14 follow the same pattern as that illustrated in Fig. 12: the soft gluon is the last (first) one to be emitted prior to (after) the shockwave. In this case too, they describe the JIMWLK evolution of one of the S-matrices from Eq. (2.10) -here, the colour quadrupole which enters in the structure of the first term there, cf. Eq. (2.14) and (2.17). Specifically, these are the four diagrams describing the evolution of the quadrupole Q (x, z, z, x) via the emission of a soft gluon which is emitted at x in the DA and reabsorbed at x in the CCA. One can check (e.g. by comparing to Eq. (A.5) in Ref. [78]) that the various factors in Eq. (5.6) are precisely as needed in order to describe this particular evolution 9 . Of course, the JIMWLK evolution of the quadrupole also involves other diagrams, in which the soft gluon is emitted/absorbed by some other pair of external legs [71][72][73]78]. These additional diagrams will be generated by the remaining contributions to the qgg cross-section, whose soft limit will be explored in the next two subsections.

Direct contributions involving the 3-gluon vertex
The discussion of the soft gluon limit for the second contribution to the qgg cross-section, Eq. (3.46), is somewhat simpler, due to the manifest symmetry of this contribution under the exchange of the two final gluons. It is therefore enough to consider only one of the two eikonal limits aforementioned, say ξ → 0. In this limit, the triple-gluon vertex in Eq. (3.24) can be approximated as It is then easily to see, by using the expressions (3.23) for the relevant effective vertex and (3.47) for the kernel, that the latter admits the same limit for all the four terms within the square brackets in Eq. (3.46), that is, This expression shows the expected soft singularity: a single pole at ξ = 0. As before, this degeneracy of the four kernels implies that only four topologies survive in the final result (among the 16 which were originally present in the r.h.s. of Eq. (3.46)): those in which the shockwave is inserted in either the final qgg state, or in the intermediate qg state. The surviving configurations are shown in Fig. 16. At large N c , each of the S-matrices associated with these four graphs is in turn built with two colour structures, corresponding to permutations of the two gluons, as shown e.g. in Eq. (3.50). Still for ξ → 0, one can also identify z = y and z = y in the arguments of the surviving S-matrices.
We are thus led to the following result (we use the conventions in Fig. 10, i.e. the k + 2 = ξq + and k + 3 = (ϑ − ξ)q + ): dσ qA→qgg+X − Q (x, z, y, x) S (z, y) + Q (x, y, y, x) S (y , y) , where Z = z − y, R = x − y and similarly for the coordinates with a bar. Both the first term within the first square bracket and the first term within the second square bracket are generated by the first graph in Fig. 16, via its decomposition at large N c , which is graphically illustrated in Fig. 17. This figure also makes clear that the four terms within the first square bracket describe the emission of a soft gluon from the quadrupole Q (x, y, y, x), such that the soft gluon is emitted at y in the DA and at y in the CCA. Similarly, the four terms within the second square bracket describe the emission of a soft gluon from the dipole S (y , y). The result (5.9) takes a more familiar form in the case when the soft gluon k 2 is not measured, meaning that one can integrate over its kinematics (which in particular allows us to identify z = z).
− Q (x, z, y, x) S (z, y) + Q (x, y, y, x) S (y , y) , This result is recognised as a piece of the JIMWLK evolution of the first term S qgqg (x, y, x, y) Q(x, y, y, x) S(y, y) in the expression (2.17) or the LO cross-section. More precisely, we have the real part of the BK equation for the dipole S(y, y) together with the relevant part of the real part of the JIMWLK equation for the quadrupole Q(x, y, y, x). To complete the (real terms in the) evolution of the latter, we also need the interference contribution in Eq. (3.52). This will be discussed in the next section. Figure 17. A redrawing of the first graph in Fig. 16 valid at large N c which exhibits the 2 colour structures hidden in this graph: each gluon is replaced by a double line denoting a quark-antiquark pair of zero transverse size. These two large-N c graphs correspond to the first terms in the first and, respectively, second squared bracket in Eq. (5.9).

Interference graphs
We finally consider the eikonal limits for the contributions to the qgg cross-section associated with interference graphs, cf. Eq. (3.52) and Fig. 11. The relevant limits for the product of effective vertices in the numerator of Eq. (3.53) read as follows: When ϑ − ξ → 0, all the four kernels in Eq. (3.52) take the same form, which reads Once again, this means that only four topologies survive in the cross-section in that limit -those in which at least one of the final gluons interacts with the shockwave and which are illustrated in Fig. 18. Still for ϑ − ξ → 0, we can also approximate the transverse coordinates as y x and y z.
Concerning the limit ξ → 0, we need to distinguish between two cases, depending upon the fact that X = x − z vanishes, or not: Clearly, the soft singularity at ξ = 0 is present only in the case where X = 0. i.e. for the two terms within the square brackets in Eq. (3.52) where the coordinates x and z get replaced by y. The eight S-matrices within the structure of these two terms can now be simply added to each other (since they are all multiplied by the same kernel) and it is easy to see that only four of them survive in the final result: those illustrated in Fig. 19. Still for ξ → 0, one can use w y (1 − ϑ)x + ϑz and y z .  To show explicit results, we shall use the momenta assignments in Fig. 11, which in particular imply k + 2 = ξq + and k + 3 = (ϑ − ξ)q + . For the case ϑ − ξ → 0, one easily finds dσ qA→qgg+X where X = x − z , Z = z − z , and we have also used y = x and y = z to replace Y → X = x − z and R → X = x − z. These alternative notations facilitate the comparison with Eq. (5.6), which refers to the direct contribution of the amplitude in Fig. 6. Both Eq. (5.6) and Eq. (5.14) express the effect of producing one soft gluon when acting with the "production" version of the JIMWLK Hamiltonian on the quadrupole Q(x, z, z, x) in the first term in the r.h.s. of Eq. (2.17) for the LO dijet cross-section 10 . In Eq. (5.6), the soft gluon is emitted and reabsorbed by the same quark, with transverse coordinates x and x in the DA and the CCA, respectively. Eq. (5.14) is the interference term where the soft gluon is emitted by the quark x in the DA and absorbed by the "other quark", at z, in the CCA. (This "other" quark is truly a component in the large-N c decomposition of the harder gluon with energy fraction ϑ.) Notice the difference in sign between Eqs. (5.6) and (5.14). Except for the sign, the respective numerical factors are identical at large N c , as they should. Besides Eqs. (5.6) and (5.14), there are two more contributions to the JIMWLK evolution of the quadrupole (here, in the sense of soft gluon production) [78]. One of them was already shown Eq. (5.9), although with somewhat different notations for the transverse coordinates 11 . The other one is the second interference term, in which one exchanges the topologies between the DA and the CCA; this can be easily inferred from (5.14) via hermitian conjugation. When the soft gluon is not measured, the sum of these four contributions reproduces the "real" terms in the B-JIMWLK evolution of the quadrupole [71][72][73].
Consider now the other eikonal limit, namely ξ → 0. Using the second line in Eq. (5.13) together with the relevant terms from Eq. (3.52), one finds dσ qA→qgg+X where Y = y − z, X = x − z , Z = z − z , and R = x − y and we have replaced w → y and z → y, as appropriate in this limit. The S-matrices within the square brackets have the right structure to describe the production of a soft gluon from the dipole S (y, y) (see Eq. (A.2) in Ref. [78]). If the soft gluon is not measured in the final state, one must integrate Eq. (5.15) over k + 2 = ξq + and k 2 ; this gives (with 2α s C F πᾱ s ) where the change of global sign occurred because Z = −(y − z). This is indeed the right sign to describe the BK evolution of the dipole S(w, y) from the product S(w, y) S(y, x) which appears with a negative sign in the third term in the r.h.s. of Eq. (2.17). (Recall that w y in the soft limit at hand.) To summarise, the discussion of the eikonal limits of our results for the qgg trijet cross-section has allowed us to reconstruct the "real" terms in the B-JIMWLK evolution of the LO cross-section (2.17) for the qg dijet. Specifically, we have identified the "real" part of the evolution equations for the first term (the product Q(x, y, y, x) S(y, y) between a quadrupole and a dipole) and for the last term (the dipole S (w, w)) in Eq. (2.17), and also for the dipoles S(y, w) and S(w, y) which appears in the two intermediate terms there. On the other hand, we have not yet encountered the evolution of the two other dipoles which enter these intermediate terms, namely S(x, y) and S(y, x). There is a good reason for that: these two other dipoles are built with quark legs which "live" on the same side of the cut, that is, they both "live" in the DA for S(x, y), and respectively in the CCA for S(y, x). Accordingly, the BK evolution of these dipoles involves graphs like those shown in Fig. 20, which represent virtual corrections to the qg dijet cross-section. The contributions of such graphs will be considered in a companion paper devoted to virtual corrections.

Collinear divergences: recovering the DGLAP evolution
Besides the soft divergences that we have just discussed, the NLO corrections to the dijet cross-section are expected to also contain collinear divergences. In momentum space, they correspond to the limit where the daughter partons emerging from a splitting have a very small relative transverse momentum (i.e. they make a very small angle θ → 0). In the transverse coordinate space, this is the limit where the transverse separation between the two daughter partons becomes arbitrarily large.
Before turning to explicit calculations, let us first recall some general features about the origin and the treatment of the collinear divergences in the framework of the hybrid factorisation (see e.g. [21]): (i) The collinear divergences are associated with (unmeasured) partons which are emitted long before the hard process, or longtime after it. By the "hard process", we mean the ensemble of the other interactions, that is, the scattering off the nuclear target (the shockwave) and the emission of additional, relatively hard (non-collinear), partons.
(ii) The "collinear" parton may scatter with the shockwave, but its scattering does not matter for the inclusive dijet cross-section at NLO.
(iii) If the collinear emission occurs prior to the hard process -in the present case, that should be a gluon emission by the incoming quark -, then the associated divergence can be reabsorbed into the DGLAP evolution of the parton distributions in the proton projectile (here, the quark distribution q f (x p , µ 2 ), as visible in Eqs. (2.16) and (4.6)).
(iv) If the collinear emission originates from one of the final partons -in this case, that could be either a gluon emission from the final quark, or the splitting of the final gluon into two gluons, or into a quark-antiquark pair -, then its treatment depends upon our definition of the measurement process. If we measure two hadrons in the final state, then the collinear divergences should be absorbed into the DGLAP evolution of the fragmentation functions for the quark, or the gluon, emerging from the hard process. If instead we measure two jets, whose definition involves an opening angle R, then any emission at small angles θ < R must be viewed as a part of the final-state jets, whereas emissions at larger angles θ > R are interpreted as NLO corrections to the hard process. In this argument, the jet angle R effectively acts as a collinear cutoff for the NLO corrections. In this remaining part of this section, we shall assume a hadronic description for the final state, for definiteness.
(v) The collinear divergences can be unambiguously separated from the infrared (or "soft") divergences associated with soft gluons. A given graph can develop both soft and collinear divergences, but the overlapping divergences cancel -as a consequence of probability conservation -when adding together real and virtual corrections. This makes it possible to disentangle soft from collinear divergences in practice and to ascribe them to the B-JIMWLK evolution of the hard process and to the DGLAP evolutions of the initial and final states, respectively.

DGLAP evolution for the initial quark
Our first example refers the NLO contribution denoted with a subscript (1) in Eq. (4.16). We recall that this contribution is obtained from the trijet cross-section (3.39) by integrating out a gluon with momentum k 3 , which can be any of the two gluons in the final state, cf. Fig. 9.
As already mentioned, the collinear regime corresponds to the situation where the transverse separation between the two daughter partons can be arbitrarily large. For the diagram in Fig. 9, there are two such limits, one for each of the two gluon emissions, as illustrated in Fig. 21. In the first case, depicted in the l.h.s. of Fig. 21, the first gluon emission is the collinear one; one then has |Y | |X |, where we recall that X = x − z and Y = y − z. This covers the case of the DGLAP evolution of the incoming quark distribution, to be discussed in this section. In the other case (the r.h.s. of Fig. 21), the opposite inequality holds: |X | |Y |. This is the case of the final state DGLAP evolution, to be discussed in the next subsection.  Consider the case depicted in the l.h.s. of Fig. 21: the first gluon emission, with longitudinal fraction ξ, is not measured (z = z) and is in the collinear regime. Using Y 2 (X ) 2 and similarly Y 2 (X ) 2 , one finds that the kernel (3.40) simplifies to (recall Eq. (3.41)) where ξ ≡ (ϑ − ξ)/(1 − ξ) denotes the splitting fraction for the second gluon emission and is the DGLAP splitting function for the process q → qg. The approximation in Eq. (6.1) holds so long as together with a similar condition involving Y and X . In writing Eq. (6.1), we have ignored the instantaneous pieces of the effective vertices which enter the numerator of K 1 (that is, we have used Eq. (3.41)). On physical grounds, it is quite obvious that these pieces cannot yield collinear singularities: they correspond to effective graphs in which the two emissions occur simultaneously. We shall shortly give a more mathematical argument in that sense. Eq. (6.1) shows that, in the collinear limit 12 , the two gluon emissions factorise from each other at the level of the kernel K 1 . Let us now check that a similar factorisation also holds at the level of the colour structure, i.e. for the S-matrices.
Since all the four terms within the squared brackets in Eq. (3.39) are now multiplied by a same kernel, it follows that only four topologies (in terms of shockwave insertions) survive in this limit. These are the topologies shown in Fig. 22, which do not involve initial-state interactions, as expected: the "collinear" gluon is emitted before the hard process, both in the DA and in the CCA. This gluon crosses the shockwave in all these graphs, but the Wilson lines describing its scattering mutually cancel by unitarity: V (z)V † (z) = 1 when z = z. Because of this cancellation, the colour structure of the four graphs in Fig. 22 is precisely the same as for the LO qg final state in Eq. (2.10).
Using Eqs. (3.39) and (6.1), one deduces the following approximation for this particular dijet cross-section in the collinear limit for the unmeasured gluon: For more clarity, we have used the same notations for the momenta of the two produced partons as in Eq. (2.10): (p + , p) and (k + , k) refer to the produced quark and gluon, respectively. Furthermore, x 1 = p + /Q + and x 2 = k + /Q + , with Q + the longitudinal momentum of the incoming proton. The original longitudinal fractions ϑ and ξ should now be evaluated as The transverse coordinate y of the intermediate quark is related to the coordinates x and z of the produced partons via Eq. (3.16), which becomes y = x 1 x + x 2 z x 1 + x 2 (and similarly for y) . (6.6) Eq. (6.4) exhibits the expected factorisation of the (unmeasured) collinear emission from the LO dijet cross-section (2.10) -here initiated by a quark with original longitudinal momentum q + (1 − ξ). The integral over z in the last line, which physically represents the integral over the transverse phase space for the collinear emission, is logarithmically divergent at large |z|. To exhibit this singularity, let us introduce a low-momentum cutoff Λ on the transverse momentum of the unmeasured gluon, meaning an upper cutoff ∼ 1/Λ on the transverse separations |y − z| and |y − z|. With this regularisation, the integral over z in Eq. (6.4) can be estimated as (below, r ≡ |y − y|) In the last equality, µ is a generic scale obeying Λ µ < 1/r. In what follows, it will be used as a renormalisation scale to subtract the collinear divergence at Λ → 0 from the genuine NLO correction. The divergence will be then absorbed into the DGLAP evolution of the quark distribution for the incoming quark.
To that aim, one should recall that the physical cross-section also includes a convolution with q f (x p , µ 2 ), as visible in Eq. (2.16) at LO and in Eq. (4.6) for a particular "real" NLO correction. For the NLO contribution at hand, the relevant convolution reads where the dots within the integrand stay for the partonic cross-section in Eq. (6.4). In the integral in the r.h.s. we changed variable x p → ξ according to (6.5); the step function comes from the condition x p ≤ 1. When using this convolution together with Eq. (6.4), one should also notice that q + (1 − ξ) = p + + k + = (x 1 + x 2 )Q + is the total longitudinal momentum of the produced dijet and thus it is independent of ξ. Hence, the whole ξ dependence of the final result is encoded in the following integral where it is understood that x ≡ x 1 + x 2 . As suggested by its notation, the r.h.s. in (6.9) is precisely the result of one "real" step in the DGLAP evolution of the quark distribution inside the proton, where the step consists in integrating out gluon radiation with virtualities q 2 within the range Λ 2 < q 2 < µ 2 .
The corresponding virtual graphs are expected to add the "plus" prescription to the splitting function P q→g (ξ) and thus remove the apparent divergence of (6.9) at ξ → 0. This will be checked in our subsequent paper (see also the discussion at the end of this subsection). By inspection of Eqs. (6.4)-(6.9), one can see that the original collinear divergence from the NLO correction has been transferred to the DGLAP evolution of the incoming quark distribution in the LO result. The finite reminder, as obtained by keeping only the second term, π ln(1/r 2 µ 2 ), in the r.h.s. of Eq. (6.7), represents a genuine NLO correction to the hard process. However, our above calculation for this correction is only correct to leading logarithmic accuracy w.r.t. the logarithm ln(p 2 s /µ 2 ), with p 2 s ≡ min(p 2 , k 2 ). (The typical value of r ≡ |y − y| is fixed by the Fourier transforms in (6.4) as r 2 ∼ 1/p 2 s .) The µ-dependence of the cross-section introduced by this and other similar corrections compensates the respective dependence of the quark distribution. This ensures that the result for the cross-section is independent of the arbitrary renormalisation scale µ at NLO.
In order to have a complete NLO result, which also includes the numerical coefficient under the log, one must go beyond our previous approximations in Eqs. (6.1) or (6.7). That is, one should exactly compute the contributions from graphs which contain collinear divergences, like the 4 graphs Fig. 22, by using dimensional regularisation together with a suitable subtraction scheme, like M S.
At this level, it is easy to understand why the graphs involving instantaneous propagators cannot generate collinear divergences. As clear by inspection of (6.4), such divergences occur via the integration over the Weiszäcker-Williams kernel Y m /Y 2 which describes a bremsstrahlung process in transverse space (recall e.g. Eq. (3.39)). This kernel is inherent in the graphs involving the regular part of the (quark or gluon) propagator, but it is absent from the instantaneous graphs. Said differently, the contribution from the instantaneous vertices do not involve enough powers of Y m and/or Y p to generate a divergence when integrating over z.
It is finally interesting to compare the graphs in Fig. 22, in which the first emitted gluon is collinear, to those in Fig. 12, where the same gluon is soft. In Fig. 22 there are no initial-state interactions, whereas in Fig. 12 the final-state interactions are missing. This is in agreement with the physical expectations that a collinear gluon must be emitted very early, well before the hard process, unlike a soft gluon, which must be the closest emission with respect to the collision. This physical distinction should allow us to unambiguously separate between soft and collinear divergences.
Yet, by inspection of these two figures, one sees that one particular graph contributes in both cases: this is the first graph in both Fig. 22 and Fig. 12. This graph has generated the term proportional to S (y, y) in Eq. (6.4), which is in fact an exact evaluation for that graph 13 -up to regularisation issues, of course. And indeed, if one takes the soft limit ξ → 0 in Eq. (6.4), one recovers one of the two terms proportional to S (w, w) in Eq. (5.5). (The other such a term comes the fourth graph in Fig. 12, in which the soft gluon is never intersecting the shockwave. We recall that, in writing Eq. (5.5), we have identified y = w and y = w.) This discussion makes clear that the first graph in Fig. 22 (or in Fig. 12) contributes to both the collinear (DGLAP) and the soft (B-JIMWLK) evolutions. In particular, it contains an overlapping, soft and collinear, divergence. This might seem to contradict our previous argument, that these two types of divergences can be separated from each other. In fact, there is no contradiction: such overlapping divergence are eventually cancelled after adding the virtual corrections.
We shall systematically study the virtual corrections in a subsequent paper. Here we merely show 13 Remember that this particular graph originates from the last term within the square brackets in (3.39), for which Eq. (6.1) becomes exact. Figure 23. The set of "real" and "virtual" graphs which count for the cancellation of the overlapping, soft and collinear, divergences associated with the first gluon emission by the quark.
in Fig. 23 the graphs relevant to the above discussion: the real graph that has already appeared as the first diagram Fig. 22 together with the virtual graphs which remove its collinear divergences in the limit where the first emission is both collinear and soft (ξ 1). The soft limit is indeed important for this argument: in general, the real NLO corrections and the virtual ones are weighted by quark distribution functions with different arguments: x p = (x 1 + x 2 )/(1 − ξ) for the real corrections, cf. Eq. (6.8), but x p = x 1 + x 2 for the virtual graphs, whose final state involves only two partons: the measured quark and gluon. However, when ξ 1, these two weighting factors become identical with each other and then it becomes possible to observe the cancellation of the overlapping, soft and collinear, divergences.
Let us show how this works for the three graphs in Fig. 23. Each of the two virtual graphs gives a contribution which is simply the product between the LO dijet cross-section (2.10) and a "tadpole" describing the virtual gluon emission. The effect of adding these contributions to the real term in Eq. (6.4) is to replace the emission kernel for the first emission, cf. Eq. (6.7), by the dipole kernel, which decreases much faster at large |z| then the original kernel in Eq. (6.7): as 1/z 4 instead of 1/z 2 . So, the would-be logarithmic singularity at large |z| disappears. The dipole kernel exhibits instead short-distance (ultraviolet) poles at z = y and z = y, but they are ultimately inocuous, as they cancel against other virtual corrections, not shown in Fig. 23.

DGLAP evolution for the final quark
Still for the topology illustrated in Fig. 9, we now consider the case where the second gluon emission -the one where the gluons carries a longitudinal momentum fraction ϑ − ξ in Fig. 9 -is not measured (z = z ) and is collinear. As already explained, this means that the transverse separation X = x − z between the two daughter partons emerging from the second splitting is much larger than the corresponding separation Y = y − z for the first emission: (X ) 2 (Y ) 2 . Under this assumption, the dominant contribution to the trijet cross-section in Eq. (3.39) comes from the last term there -the one obtained after the double replacement x, z → y & x, z → y -since the corresponding kernel is the only one not to be suppressed at large (X ) 2 . Not surprisingly, the graphs which survive in this limit are those which are void of final state interactions, as illustrated in Fig. 24. We therefore expect the collinear divergence generated by these graphs to express the DGLAP evolution of the final-state gluon in the leading-order dijet cross-section (2.10). Let us verify that this is indeed the case.
The kernel for the surviving terms has, clearly, the same expression as displayed in Eq. (6.1), where the longitudinal momentum fractions must now be evaluated as The main difference w.r.t. the previous subsection is that the kernel in Eq. (6.1) now applies only to the last term within the square brackets in Eq. (3.39). After also using z = z , we find that Eq. (6.4) gets replaced by dσ qA→qg+X For reasons to shortly become clear, it is convenient to consider z 1 ≡ 1 − ξ -the splitting fraction of the measured quark at the second emission vertex -as an independent variable, on the same footing as the external variables x 1 and x 2 . Then the longitudinal fractions ξ and x p and the transverse coordinate y should be understood as (recall Eqs. (3.16) and (6.11)) together with a similar expression for y.
It is furthermore useful to change two of the integration variables, from x and x to y and y (this will facilitate the comparison with the LO result in Eq. (2.10)). Recalling that z = z , one finds p · (x − x) = 1 z 1 p · (y − y), x − z = 1 z 1 (y − z ) . (6.14) When expressed in terms of these new variables, w = (1 − ξ)y + ξz becomes independent of z , so the integral over z factorizes. After also convoluting with the initial quark distribution, as shown in the l.h.s. of Eq. (6.8), and changing the respective integration variable from x p to z 1 according to Eq. (6.13), one obtains the following expression for the collinear singularity encoded in this particular NLO correction to the dijet cross-section dσ pA→qg+X (1) rNLO,2 dp + d 2 p dk + d 2 k dz 1 z 3 1 x p q f (x p , µ 2 ) 4α s C F (2π) 6 (q + ) 2 P q→g (ξ) × y, z, y, z e −ip·(y−y)/z 1 −ik·(z−z) (y − z) · (y − z) (y − z) 2 (y − z) 2 × Q(y, z, z, y) S(z, z) − S(y, z) S(z, w) − S(w, z) S(z, y) + S (w, w) where we have also used x 1 /(1 − ξ) = x p z 1 and P q→g (1 − z 1 ) = P q→q (z 1 ). Eq. (6.12) exhibits a factorised structure, as expected: it is the product of the LO dijet crosssection in Eq. (2.10) (but for final momenta p 1 = p/z 1 and k) times the probability for an unobserved emission in the final state. This probability contains a collinear singularity similar to that visible in Eq. (6.7), which here should be absorbed into the renormalisation of the fragmentation function of the final quark. Specifically, the singular piece of Eq. (6.12) can be written as dσ pA→qg+X (1) rNLO,2 dp + d 2 p dk + d 2 k coll = dz 1 z 3 1 x p q f (x p , µ 2 ) dσ pA→qg+X LO d 3 p 1 d 3 k p 1 =p/z 1 D g/q (z 1 , µ 2 ) , (6.16) where D g/q (z 1 , µ 2 ) ≡ α s C F π P q→q (z 1 ) ln µ 2 Λ 2 . (6.17) is recognised as the contribution of the first step in the DGLAP evolution of the quark-to-gluon fragmentation function. To summarise, Eqs. (6.15)-(6.17) describe the beginning of the DGLAP evolution of the quark fragmentation function in Eq. (2.18). The corresponding evolution of the gluon fragmentation function will be discussed in the next two sections.
6.3 DGLAP evolution for the final gluon: g → gg splitting We now consider the collinear limit for the second topology yielding real NLO corrections to quarkgluon production, the one denoted by the subscript (2) in Eq. (4.16) and which is illustrated in Fig. 7. In this case, a collinear divergence can be generated only by the second emission (see the discussion Figure 25. The four graphs surviving in the collinear limit for the g → gg splitting (the unmeasured gluon is depicted in blue). These graphs are void of final-state interactions.
towards the end of this section), in which case it is associated with the DGLAP evolution of the fragmentation function for the final gluon.
So let us consider the collinear limit for the second emission in Fig. 7. This is the limit where the transverse separation Z = z − z between the two final final gluons is much larger than the separation R = x−y between the final quark and the intermediate gluon: Z 2 R 2 . With reference to Eq. (3.46), this hierarchy entails two important simplification. The first one refers to the energy denominator and is fully similar to that discussed in the previous section: it implies that the only surviving graphs are the four graphs without final-state interactions, shown in Fig. 25. These graphs correspond to the last term within the square brackets in Eq. (3.46) -the only one not to be suppressed at large Z 2 . The second simplification refers to the product (3.48) of effective vertices: in Eq. (3.46), this product is multiplied by R m Z n R p Z q . When computing the graphs in Fig. 25, the integral over the coordinate z of the unmeasured gluon (z = z ) can be factorised as where the approximate equality holds for the singular piece alone: in the collinear limit at hand, the differences Z = z − z and Z = z − z can be arbitrarily large, whereas |z − z| is fixed to a value ∼ 1/k ⊥ by the Fourier transform for the measured gluon; so, in evaluating the singular piece of (6.18), one can approximate Z Z and at the same time restrict the integral over z to Z 2 > (z − z) 2 . The logarithmic divergence at Z 2 → ∞ has been regulated by a momentum cutoff Λ 2 , as in Eq. (6.7).
Eq. (6.18) shows that, for the purpose of extracting the collinear singularity, one can replace Z n Z q → (Z · Z)δ nq /2 within the integrand of Eq. (3.46). This allows us to simplify the tensorial structure of the product of effective vertices (cf. Eq. (3.48)): This is analog to Eq. (6.15) from the previous section and its subsequent discussion is also similar: the collinear divergence can be absorbed in one step of the DGLAP evolution of the gluon-to-gluon fragmentation function: D g/g (z 2 , µ 2 ) ≡ α s N c π P g→g (z 2 ) ln µ 2 Λ 2 . (6.25) To conclude this section, let us explain why the first emission cannot contribute a collinear divergence for this particular topology. The first emitted gluon is an intermediate gluon, which decays into two other gluons prior the final state. One of these daughter gluons is measured in the final state (z = z), and the other one is not (z = z ). Hence the coordinates, y and y, of the intermediate gluon are different in the DA and the CCA, respectively (as also manifest in Eq. (6.20)). Hence, the would-be collinear limit for the first emission, that is, R 2 Z 2 , cannot generate a collinear divergence in the integral over the transverse coordinate z of the unmeasured daughter gluon.
One can similarly understand that the interference graphs responsible for the piece denoted by the subscript (3) in Eq. (4.16) (see Fig. 11) do not generate collinear divergences either. This is consistent with the fact that the LO DGLAP evolution admits a probabilistic picture and could not accommodate such interference effects.

(A.9)
The multi-parton bare Fock states are obtained by acting with the relevant creation operators on the bare vacuum state. In this paper, we use both the 3-momentum representation k = (k + , k) and the mixed representation (k + , x), as obtained via the Fourier transform from transverse momenta to transverse coordinates. Let us present here a few representative examples.