Classical Gluon Production Amplitude for Nucleus-Nucleus Collisions: First Saturation Correction in the Projectile

We calculate the classical single-gluon production amplitude in nucleus-nucleus collisions including the first saturation correction in one of the nuclei (the projectile) while keeping multiple-rescattering (saturation) corrections to all orders in the other nucleus (the target). In our approximation only two nucleons interact in the projectile nucleus: the single-gluon production amplitude we calculate is order-g^3 and is leading-order in the atomic number of the projectile, while resumming all order-one saturation corrections in the target nucleus. Our result is the first step towards obtaining an analytic expression for the first projectile saturation correction to the gluon production cross section in nucleus-nucleus collisions.


I. INTRODUCTION
Complete understanding of heavy ion collisions is impossible without a clear picture of the initial stages of the collision, i.e., the processes leading to the creation of quarks and gluons which later on thermalize forming quarkgluon plasma (QGP). The distribution of quarks and gluons produced in these early-time processes, known as the initial condition for the QGP formation, is the fundamental building block of heavy ion theory. Questions concerning thermalization of the quark-gluon medium and the determination of the initial conditions for the hydrodynamic evolution of the QGP can not be answered in a fully satisfactory manner without the qualitative and quantitative knowledge of the production mechanism for the initial-state quarks and gluons in a nuclear collision.
In the framework of saturation physics (see the reviews [1][2][3][4][5][6][7] and the book [8]) the leading-order contribution to gluon production is given by the classical gluon fields [9][10][11][12][13] of the McLerran-Venugopalan (MV) model [14][15][16]. Classical gluon fields in heavy ion collisions resum powers of α 2 s A 1/3 1 and α 2 s A 1/3 2 [17], where α s is the strong coupling constant, while A 1 and A 2 are the atomic numbers of the two nuclei (henceforth referred to as the projectile and the target). Since the saturation scales squared of the two nuclei are proportional to these parameters, Q 2 s1 ∼ α 2 s A 1/3 1 and Q 2 s2 ∼ α 2 s A 1/3 2 , we can write down the quasi-classical single-gluon production cross section as where B ⊥ is the impact parameter between the two nuclei, b ⊥ is the transverse position of the produced gluon with respect to the center of the target nucleus, while k ⊥ is the transverse momentum of the produced gluon with k T = | k ⊥ |. The expansion in the powers of α 2 s A 1/3 1 and α 2 s A 1/3 2 corresponds to expansion in the powers of Q 2 s1 /k 2 T and Q 2 s2 /k 2 T , Note that due to projectile-target symmetry (resulting in Q s1 ↔ Q s2 symmetry) we have c n,m = c m,n . (The k Tdependence on the right-hand side of Eq. (1) also enters through powers of ln(k T /Λ) where Λ is the infrared (IR) cutoff of each nucleon: the powers of ln(k T /Λ) are not shown explicitly above. The coefficients c n,m from Eq. (2) are, in fact, polynomials in ln(k T /Λ).) At the moment we do not have an analytic expression for the function f (Q 2 s1 /k 2 T , Q 2 s2 /k 2 T ) in Eq. (1). This function was extensively studied numerically in [18][19][20][21][22][23][24]. Still it appears desirable to attain a better handle on the analytic form of this function. Apart from the general advantage of having an analytic solution, knowing the function f (Q 2 s1 /k 2 T , Q 2 s2 /k 2 T ) should greatly facilitate the inclusion of small-x evolution corrections [25][26][27][28][29][30][31][32][33][34] along with the running-coupling corrections [35][36][37][38][39] into the gluon production cross section. These corrections are essential for realistic phenomenological applications.
Let us briefly summarize what is known about the function f (Q 2 s1 /k 2 T , Q 2 s2 /k 2 T ). The leading-order result (the coefficient c 1,1 in Eq. (2)) was obtained in [9,10,12], reproducing the earlier results of [25,40]. The case of protonnucleus (p + A) collisions, defined as the leading-order in Q 2 s1 term in the expansion on the right-hand side of Eq. (2), was solved in [41] (see also [42,43]), yielding the coefficients c 1,n (and, due to target-projectile symmetry, c n,1 ) for any positive integer n. This is all we presently know analytically about the coefficients c n,m .
An ansatz for the full solution of the classical gluon production problem was proposed by one of the authors in [44]. A variational approach to the problem was attempted in [45]. While consistent with our knowledge of coefficients c 1,n , neither of these results can be verified further due to our lack of knowledge of the coefficients c n,m for n, m ≥ 2. In phenomenological applications one often employs the k T -factorization formula involving unintegrated gluon distributions [46][47][48]: while the gluon production cross section in p + A collisions (the lowest-order in Q 2 s1 terms in Eq. (2)) does lead to the k T -factorization formula [49,50], it is not clear whether k T -factorization holds beyond the p + A approximation, again due to our lack of knowledge of c n,m for n, m ≥ 2. Moreover, numerical simulations of the classical gluon production [18][19][20][21][22][23][24] for nucleus-nucleus (A + A) collisions appear to rule out the k T -factorization ansatz, suggesting that this factorized result is only valid for p + A.
The goal of the current project is to determine the coefficients c 2,n (for n ≥ 2). As one can see from Eq. (2), to obtain c 2,n coefficients one has to expand the gluon production cross section to the second order in Q 2 s1 keeping all orders of Q 2 s2 . This means that one has to allow two nucleons in the projectile nucleus to participate in the interaction, while allowing all nucleons in the target nucleus to interact. Formally one can think of this as working in the regime where α 2 s A correction to gluon production. Note that one still has A 1 ≫ 1 such that the O(α 3 s ) contribution to the cross section, where only one of the projectile nucleons interacts, is small due to it being suppressed by a power of A 1/3 1 . As we will detail below, having two nucleons interact in the projectile nucleus results in the gluon production cross section consisting of two contributions (see the top two diagrams in Fig. 2 below): (i) each of the two nucleons interacts both in the amplitude and in the complex conjugate amplitude (or one nucleon interacts only in the amplitude and another nucleon interacts only in the complex conjugate amplitude), or (ii) one nucleon interacts both in the amplitude and in the complex conjugate amplitude while the other nucleon interacts only in the (complex conjugate) amplitude. In the case (i) the scattering amplitude is O(g 3 ), such that the contribution to the cross section is ∼ |g 3 | 2 ∼ α 3 s . In the case (ii) the amplitude is O(g 5 ), while the complex conjugate amplitude is O(g), such that the cross section contribution is ∼ g 5 g ∼ α 3 s . (As we will argue below, the case where the amplitude is O(g 4 ) and the complex conjugate amplitude is O(g 2 ) is included by using retarded Green functions in the O(g 3 ) and O(g 5 ) amplitudes in this quasi-classical calculation.) In the present paper we calculate the scattering amplitude in case (i). While its square would give a contribution to the desired gluon production cross section at order-Q 4 s1 , the complete expression for the cross section can only be obtained if one includes the contribution of case (ii) as well. Calculation of the contribution (ii) would involve the O(g 5 ) amplitude, which appears to be significantly more involved and is left for future work. Our calculation is performed in the light-cone gauge of the projectile, with the final results given in Eqs. (36) and (40).
The paper is structured as follows. In Sec. II we set up the problem by first reproducing the lowest-order (p + A) gluon production calculation and then by outlining the main ingredients needed to complete the calculation of the coefficients c 2,n in Eq. (2). The elements of the calculation, along with the main results, are presented in Sec. III. We conclude in Sec. IV.

II. THE SETUP
Consider the high-energy scattering of a projectile nucleus on a target nucleus. Defining the light-cone variables as v ± = (v 0 ± v 3 )/ √ 2 with the x 3 -axis being the collision axis, we choose the projectile to be moving along the x + light-cone direction, and the target moving along the x − direction. All of our calculation will be performed in A + = 0 light-cone gauge. Transverse plane vectors are denoted as v ⊥ = (v 1 , v 2 ), such that the full 4-vector is v µ = (v + , v − , v ⊥ ) in the (+, −, ⊥) notation.
A. Classical gluon production in the p + A approximation To introduce our formalism, let us begin by outlining the calculation of the single gluon production cross section in p + A collisions in the quasi-classical approximation. Working in the approximation where α 2 s A is required we only need to include diagrams where one nucleon from the projectile interacts with the target. The diagrams contributing to the gluon production cross section in p + A are shown in Fig. 1. The target nucleus moving along the x − -axis generates the shock wave, shown in this paper as a vertical band. Diagrammatically this vertical band represents multiple interactions with the field of the shock wave, which happen over a very short time interval around x + = 0. Since only one nucleon in the projectile is involved in the interaction, we model it with a single quark, shown in Fig. 1, by a horizontal solid straight line. The spectator quarks are not shown explicitly in Fig. 1 (along with the spectator nucleons if the projectile is a nucleus): below, for other processes, we also do not show the spectators explicitly. The produced gluon can be emitted either before the interaction with the shock wave (left diagram in Fig. 1) or after the interaction (right diagram in Fig. 1). Emission during the passage of the projectile through the shock wave is suppressed by a power of center-of-mass energy. Gluon emission from within the shock wave is also suppressed in the A + = 0 light-cone gauge we are using. Normalizing the incoming quark in the projectile as a Wilson line, we write the following expression for the scattering amplitude resulting from the diagrams in Fig. 1: The amplitude in Eq. (3) includes a Fourier-transform from transverse momentum space to transverse coordinate space. It also involves the adjoint and fundamental Wilson lines describing respectively the propagation of a high-energy gluon and quark through the target gluon field A µ . Here T a and t a are correspondingly the adjoint and fundamental SU(N c ) generators. In arriving at Eq. (3) we used the Fierz identity U ab V t b = t a V when evaluating the right diagram in Fig. 1 in order to put the contributions of both diagrams into a similar form.
Performing the Fourier transform in Eq. (3) we obtain The gluon production cross section is given by (see e.g. [8]) where the summation over colors and polarizations of the final-state particles along with the averaging over the quantum numbers of the initial-state particles are implicitly implied. The angle brackets . . . denote averaging in the target nucleus wave function.
Substituting the amplitude Eq. (6) into Eq. (7) yields [41] dσ where the gluon dipole S-matrix is defined by In the quasi-classical MV/Glauber-Mueller (GM) approximation it is given by [51] is the square of the gluon saturation scale with T ( b ⊥ ) the nuclear profile function and Λ the IR cutoff of each individual nucleon (Λ ∼ Λ QCD ). Eq. (10) resums all the multiple rescatterings in the target nucleus: thus, in the MV model, it resums all-order saturation corrections in the target.
Note a particular convenience of the form of the amplitude given in Eq. (6), with the fundamental Wilson line V placed to the left of the fundamental SU(N c ) generator in both terms there: when we square the amplitude, the fundamental Wilson line V is multiplied by its hermitean conjugate V † giving an identity. Below when calculating diagrams we will always cast them in the same form, which would make all fundamental Wilson lines vanish when we square the amplitude.

B. First Saturation Correction in the Projectile: Diagram Types and Calculational Simplifications
The goal of this project is to calculate the first projectile saturation correction to the cross section in Eq. (8). That is, we need to find the order-α 2 s A 1/3 1 correction to (8): this is the leading in A 1 part of the order-α 2 s correction. To get the leading-A 1 contribution, the correction must include an interaction with another nucleon in the projectile nucleus. Hence we see that we need to find order-α 2 s correction to (8) involving one other nucleon in the projectile. The main types of the diagrams we need to calculate are shown in Fig. 2. There by two horizontal solid straight lines we show two quarks from two different nucleons in the projectile nucleus, again suppressing the spectators. The diagrams in Fig. 2 represent the amplitude squared which contributes to the cross section: the solid vertical line denotes the final-state cut. The cross labels the measured gluon.
The three types of projectile saturation corrections to gluon production shown in Fig. 2 are (i) the square of the order-g 3 amplitude; (ii) the interference between the order-g 5 amplitude and the leading-order (order-g) amplitude from Eq. (6); and (iii) the interference between the order-g 4 amplitude and the order-g 2 amplitude. (In this power counting we are assuming that the interaction with the shock wave is order-1, that is α 2 s A 1/3 2 ∼ 1, and are counting only powers of g arising from the vertices shown explicitly in Fig. 2.) The order-α 3 s diagrams in which one of the nucleons is a spectator and does not participate in the interactions are suppressed by A 1/3 1 and are neglected in our analysis. Note that diagram (iii) has two gluons in the final state: since we are calculating the single-inclusive production cross section, having more than the measured gluon in the final state is allowed.
Below we calculate the order-g 3 amplitude that enters diagram (i) in Fig. 2, thus constructing the order-α s correction to the amplitude in Eq. (6) which is enhanced by the leading power of A 1/3 1 . The order-g 5 amplitude from diagram (ii) is left for future work.
The situation with the diagram (iii) is more subtle. First let us note that, since we are working in the classical MV model, the diagrams we consider for gluon production also correspond to the diagrams contributing to the classical gluon field with the two colliding nuclei providing the source current [9,10,12,[14][15][16]. Hence one can think of these diagrams as graphically representing the solution of the Yang-Mills equations [12,17]. In such case each Feynman propagator would be replaced by the retarded Green function. While such a replacement is relatively straightforward for gluon propagators, it is less clear what this prescription means for quark propagators in the time-ordered picture we employ here: for instance, the solution of Yang-Mills equations can only have the produced gluon emitted by the projectile quark after the interaction with the shock wave (see e.g. [12]), whereas a naive application of perturbation theory allows one to draw a number of a priori non-zero graphs with the produced gluon emitted before the shock wave interaction (see Fig. 7 below).
The diagrams representing different types of saturation corrections in the projectile nucleus to the gluon production cross section.
To clarify this issue we need to perform a detailed diagrammatic analysis. First let us review how the gluon Feynman propagators become retarded Green functions. As an example consider the top two diagrams shown in Fig. 3 where regular Feynman propagators are implied for uncut lines. The difference between the top two diagrams in Fig. 3 is due to one of the vertices involving the gluon carrying momentum l being moved across the cut in the right diagram. The rest of the diagram is the same in both cases. Concentrating on the propagator of the gluon carrying momentum l and suppressing the rest of the diagrams' contribution we see that adding the two graphs gives [52] −i D µν (l) 3. An example of diagrams which add up to convert a gluon Feynman propagator into a retarded Green function. The bold arrow on the gluon line in the last diagram denotes a retarded gluon Green function, with the arrow's direction pointing in the direction of light-cone time flow, that is, indicating that x + 2 > x + 1 . where is the numerator of the light-cone gauge gluon propagator with η µ = (0, 1, 0 ⊥ ) in the (+, −, ⊥) notation, such that A + = η · A = 0 is the gauge condition. For an on-shell (l 2 = 0) gluon, the numerator of the gluon propagator in Eq. (12) can be written as a sum over physical (transverse) gluon polarizations with the polarization 4-vector in the light-cone gauge and ǫ λ ⊥ = −(1/ √ 2)(λ, i). (Note that the last term on the right hand side of (12) vanishes if l 2 = 0.) In arriving at Eq. (11) we have also used the fact that the quark-gluon vertex changes sign when carried across the cut due to complex conjugation.
From Eq. (11) we conclude that contributions of top two diagrams in Fig. 3 can be found by only calculating the left diagram with the retarded gluon Green function: this conclusion is illustrated in the bottom diagram in Fig. 3, with the retarded gluon Green function denoted by a bold arrow on the gluon line. Note also that the retarded Green function arising in Eq. (11) implies that x + 2 > x + 1 in the notation shown in Fig. 3, such that the gluon is first emitted by the quark line, and then it enters the triple gluon vertex leading to the production of the measured gluon. The arrow in the last diagram of Fig. 3 indicates this direction of time flow.
So far we have considered an example of just two diagrams which combine to give us a retarded gluon Green function. However, the statement that by moving an end of a gluon line across the cut and by adding that contribution to the original diagram we would obtain a retarded Green function for that gluon is valid in general, at this classical level. To demonstrate this explicitly we need to consider many different types of diagrams: this is done in Appendix A. Note also that our ability to use a retarded gluon Green function in calculating the amplitude should not depend on what goes on in the complex conjugate amplitude: while naively our argument in Fig. 3 seems to depend on the absence of final-state (post-shock wave) interactions in the complex conjugate amplitude, the argument is in fact true in general, as long as we are working at the classical level, that is, calculating order-α 3 s cross-section contribution involving two interacting projectile nucleons (quarks). We refer the reader to Appendix A for details.
The diagram in the top left panel of Fig. 3 is in the (i)-class by the classification presented in Fig. 2, while the diagram in the top right panel of Fig. 3 is in the (iii)-class. We see that the contribution of the type-(iii) diagram is included in the type-(i) diagram by using a retarded gluon propagator in the latter. Again this conclusion is true in general: type-(iii) diagrams are included in the calculation of type-(i) and type-(ii) diagrams if we use retarded gluon propagators. (3) FIG. 4. Sample diagrams with a gluon exchange between the projectile quarks. Now let us consider diagrams where the projectile quarks exchange a gluon with each other (the gluon may or may not go through the shock wave). An example of several such diagrams is given in Fig. 4, where graphs labeled (1), (2) and (3) are different from each other only by the placement of the gluon exchanged between the projectile quarks. Once again, the arrow on this gluon line denotes a retarded gluon Green function, with the direction of the arrow indicating the direction of the resulting (light-cone) time ordering. We see that each diagram in Fig. 4 represents a sum of two diagrams (akin to Fig. 3): each graph of Fig. 4 implies a sum of itself (with the Feynman gluon propagator for the gluon with the arrow on it) and the contribution where the lower quark-gluon vertex of the same gluon is carried across the cut (to the region before the shock wave interaction in the complex conjugate amplitude), resulting in the retarded Green function for the gluon marked by an arrow. Note that the time-ordering of the retarded Green function for the gluon in the graph of Fig. 4 again points towards the produced gluon, similar to the case illustrated in Fig. 3. To analyze the diagrams in Fig. 4 we need the following additional observation: moving a retarded gluon Green function across the cut flips the sign of the contribution. This is illustrated by a simple example in Fig. 5 (see Appendix A for more details), where the cancellation is valid independent of the interactions which may happen after the shock wave interaction (i.e., for x + > 0 shown by the shaded region in Fig. 5) on both sides of the cut.
Employing the cancellation of Fig. 4, we see that diagram (1) in Fig. 4 is canceled by the contribution to diagram (3) coming from the x + 2 > x + 1 region, Due to the difference in color factors, such cancellation does not happen between the diagram (2) in Fig. 4 and the diagram (3) with x + 2 < x + 1 : instead, if we write the contribution of diagram (2) factoring out the color factor as (2) = t a t b M , we get This conclusion is illustrated in the lower right diagram of Fig. 4: the sum of diagrams (1), (2) and (3) is simply diagram (2) with the commutator instead of the product of the fundamental generators in the amplitude. Once again, while we have illustrated the point by a simple example, it is valid in all cases: contributions of diagrams with the single or double gluon exchange between the two valence quarks can all be found by simply calculating a small subset of those graphs replacing products of t a 's in them by commutators. The details of how this happens and which diagrams survive in the end are given in Appendix A. At order-g 3 the surviving diagrams are B 2 , B 9 and B 11 from Fig. 7 below (indeed diagram B 2 is the amplitude to the left of the cut in the graph (2) of Fig. 4), which have to be calculated with the commutators instead of the products of fundamental generators. As we will see below, B 2 = 0 (even with the commutator), such that only produced gluon emissions by projectile quarks after the shock wave interaction, that is for x + > 0, survive. This observation completes the analogy between the Feynman diagrams in the shock wave formalism and the diagrams representing the classical gluon fields in the MV model. We would like to stress that to establish this diagrams versus classical fields analogy, and simplify our calculations in the process, we had to utilize the fact that the amplitude we are calculating would have to be squared to obtain the cross section. Hence we have absorbed some contributions from the complex conjugate amplitude into our amplitude by employing the retarded Green functions and, for some graphs, commutators. Therefore, strictly-speaking, below we will not be calculating a standard Feynman-diagram amplitude, but a somewhat modified amplitude with retarded "propagators" and commutators: the square of this amplitude would still give the gluon production cross section.
Let us summarize the main conclusions of this Section. First of all, to find the first saturation correction in the projectile, we only need to calculate O(g 3 ) and O(g 5 ) gluon production amplitudes with the retarded gluon Green functions. The time-ordering of the retarded "propagators" should be such that the gluon would be emitted first, and then would participate in an interaction ultimately leading to the production of the tagged gluon. Secondly, the contributions of the diagrams with the gluon exchanges between the projectile nucleons (quarks) can be more economically constructed by calculating a subset of those graph with the commutators of fundamental color matrices.

A. Graphs A, B, and C
The diagrams contributing to the gluon production amplitude at O(g 3 ) that include interactions with both nucleons are shown in Figs. 6 and 7. The diagrams from Fig. 6 are labeled A i with i = 1, . . . , 7 and will be referred to as the A-graphs. Similarly the graphs in Fig. 7 are labeled B i with i = 1, . . . , 12 and will be called the B-graphs. In addition to all the B-graphs one has to consider the diagrams similar to those in Fig. 7 but with the two projectile quarks interchanged: these will be labeled C i with i = 1, . . . , 12 and referred to as C-graphs. To calculate the diagrams in Figs. 6 and 7 we will follow the rules established in the previous Section. All gluon propagators are retarded, such that time flows towards the measured gluon in the diagram. This means that the retarded propagator for the gluon exchanged between the two quarks in Fig. 7 is such that the lower quark-gluon vertex comes earlier than the upper one. The projectile quark lines are normalizes as Wilson lines (that is, each line is divided by 2P + where P + is the large "+" momentum in the line).
Our ultimate goal is to obtain an amplitude in coordinate space, in a form similar to Eq. (6). The reason for that is the relative ease with which the correlators of Wilson lines are calculated in the transverse coordinate space, particularly in the MV model. For the diagrams in Figs. 6 and 7 we would have to perform Fourier transforms into transverse coordinate space, similar to the transition from Eq. (3) to Eq. (6) above. Our standard notation will be to label the projectile quarks 1 and 2, such that their transverse positions are b 1⊥ and b 2⊥ . Note also that a priori the two quarks in the x + -direction moving projectile have different x − positions, labeled here as b − 1 and b − 2 (see diagram A 1 in Fig. 6). Since the projectile is also a large nucleus, the difference b − 1 − b − 2 , while suppressed by a power of energy, is enhanced by a power of the projectile atomic number A 1 [11,13,17]: in the following we will keep the difference b − 1 − b − 2 non-zero throughout the calculation while remembering that it is sub-eikonally small, and will put it to zero at the end of the calculation where allowed. The difference b − 1 − b − 2 will serve as a regulator in the longitudinal Fourier transform. This is in accordance to the standard regularization of singularities in the MV model [11,13,17].
To obtain an amplitude dependent on b − 1 and b − 2 we need to perform a longitudinal Fourier transform integrating over the "+" component of the momentum exchanged between the projectile quarks in the diagrams of Figs. 6 and 7. Since we are working in the A + = 0 light-cone gauge with the numerator of the gluon propagator given in Eq. (12), we have to specify the regularization of the light-cone poles at l + = 0 in order to perform this longitudinal Fourier transform. Each way of regulating the light-cone pole is equivalent to specifying a particular sub-gauge within the The connected diagrams without a triple-gluon vertex. The remaining twelve graphs of this type are obtained by swapping the Wilson lines in the diagrams shown, Ci(b1, b2; 1, 2) = Bi(b2, b1; 2, 1). All gluon propagators are retarded, with the time flowing in the direction of the produced gluon.
light-cone gauge. In our calculation we will use the principal value (PV) prescription, where and The PV prescription corresponds to the field of the projectile nucleus moving in the x + direction satisfying the sub-gauge condition (see e.g. [53]). Other possible sub-gauge choices include requiring that the field of the projectile nucleus obeys the A ⊥ (x − → +∞) = 0 condition [52,53]. This corresponds to [17] if the momentum l flows from the index µ to the index ν along the gluon line. The corresponding polarization 4-vector is One may also employ the A ⊥ (x − → −∞) = 0 sub-gauge condition, which results in reversing the sign of iǫ in Eqs. (18) and (19). The reason we are going to use the PV sub-gauge is explained in detail in Appendix B. It turns out that while the A ⊥ (x − → +∞) = 0 and A ⊥ (x − → −∞) = 0 sub-gauges are very useful for many classical gluon field calculations [11,17,41,54,55], in our shock-wave target setup the PV prescription is more economical. Namely, as detailed in Appendix B, using the A ⊥ (x − → +∞) = 0 or A ⊥ (x − → −∞) = 0 gauge choices leads to the need to include new diagrams, which are not included in Figs. 6 and 7 (and are not part of the C-graphs). See for instance the first two diagrams in Fig. 18 of Appendix B. These diagrams were considered before in [41]. In the shock-wave formalism these diagrams are probably classified as higher-order corrections to the interactions with the shock wave: namely, if one takes the leading-order gluon production diagrams from Fig. 1, and considers an order-α s correction to the interactions of either the quark or the gluon with the shock wave, one would then obtain an order-g 3 contribution to the gluon production amplitude. Indeed such corrections are not going to be enhanced by an extra A 1/3 1 (in the cross section): however, shock wave interaction corrections to the lowest-order graphs from Fig. 1 involving a quark from another projectile nucleon would give an order-g 3 contribution enhanced by a power of A 1 , that is, they would be comparable to the A, B and C graphs. In Appendix B we show that while in the A ⊥ (x − → +∞) = 0 and A ⊥ (x − → −∞) = 0 sub-gauges such diagrams are important, most of their contributions are zero in the PV sub-gauge. The remaining shock wave corrections contributions which are non-zero in the PV sub-gauge cancel in the amplitude squared, and, hence, can also be neglected in the calculation. We will therefore proceed by calculating the Fourier transform for all diagrams using the PV sub-gauge of the light-cone gauge.
We performed all diagram calculations treating the projectile quark lines both as the regular high-energy quarks and as the light-cone Wilson lines from Eq. (5) (only with non-zero b − ). In the PV sub-gauge both results are identically equal. However, in the A ⊥ (x − → +∞) = 0 and A ⊥ (x − → −∞) = 0 sub-gauges we found that the Wilson-line approximation does not give the right answer for the B and C graphs. For more details we refer the reader to Appendix B.
Since presenting a detailed calculation of each diagram is rather tedious, and would make the paper difficult to read, we will first work out one sample diagram, and then state the answers for the rest of the graphs. It appears that the diagram whose calculation illustrates most of the issues relevant to computing A, B and C graphs is A 2 . The diagrams is shown with all the momentum, coordinate, polarization, Lorentz-index and color labels in Fig. 8.
Using Fig. 8 and treating the projectile quarks as light-cone Wilson lines we write the contribution of this diagram as The integrals over x + 1 and x + 2 , coming from the projectile quarks Wilson lines, along with the integral over y + , are regulated by the exponentials like e ±ǫ x + to make them convergent at their respective infinities. This regularization is consistent with the regularization of Feynman propagators. Note that in Eq. (20) we are first Fourier-transforming the diagram into coordinate space, then integrating over x + 1 , x + 2 and y + with the proper light-cone time ordering and over x 1⊥ and x 2⊥ . The shock wave is at x + = 0. Due to the nature of the Fourier transform, the "+" components of momentum are conserved in all interaction vertices, including the interactions with the shock wave. Conversely, the "−" momentum component is not conserved at any of the vertices, such that e.g. (k − l) − = k − − l − , just like in light-cone perturbation theory (LCPT) [56], though our internal lines are not on mass shell. The outgoing gluon brings in a Fourier factor of e i k − y + , where k − = k 2 ⊥ /(2k + ) since it is on mass shell. Note that k + > 0. The gluon interaction with the shock wave brings in a factor of (2k + ) U ab x ⊥ g µν for a gluon line with light-cone momentum k + and transverse coordinate x ⊥ , according to the standard rules of the eikonal approximation [27]. Similarly, eikonal interaction of a quark with a shock wave yields (2P + ) V b ⊥ for a quark line with light-cone momentum P + and transverse coordinate b ⊥ : as we have mentioned before, the factor of (2P + ) is removed by our normalization of the external quark lines. To simplify the calculations, we have replaced the numerators of the gluon propagators by polarization sums in Eq. (20): this is justified since (cf. Eq. (12)) where the contribution of the "instantaneous" last term on the right vanishes in A 2 . Note that we are using retarded Green function regularization of all gluon propagators, in agreement with the discussion above and the arguments presented in Appendix A.
Integrating Eq. (20) over x + 1 , x + 2 and y + and employing Eq. (17) yields Since the expression in the square brackets of Eq. (22) is independent of the "−" components of momenta, we can where we have also simplified the expression in the square brackets using Eq. (17) again and defined α = l + /k + . Integration over α in Eq. (23), while a little lengthy, can be straightforwardly performed using residues. In the spirit of regulating the α-integral with b − 2 − b − 1 , we will use b − 2 − b − 1 to tell us which way to close the contour for divergent terms in the α-integral, and drop it in the exponent afterwards. For the convergent terms in the α-integral we can put b − 2 − b − 1 to zero from the start. To be more specific, consider the following integral (not present in Eq. (23)): First of all, without the exponential containing b − 2 −b − 1 the integral is divergent, and it is regulated by the exponential. In the last step in Eq. (24) the Sign function, so we leave it as is. The above is the standard procedure when calculating observables in the MV model [11,13,17]. After squaring the amplitude we would have to average the resulting cross section in the projectile nucleus wave function, which would include integrating over all b − 1 and b − 2 . This is left for future work since, as explained above, we are not ready to calculate the gluon production cross section having constructed only the O(g 3 ) amplitude in this paper.
After carrying out the α-integral Eq. (23) reduces to where we have defined Finally noticing that we can further simplify Eq. (25) to the form given below in Eq. (28b).

Results
Before performing the transverse Fourier transforms, let us first list the expressions for the contributions of all the diagrams in Figs. 6 and 7 in transverse momentum space. All the diagram values below are given in the notation where the fundamental Wilson lines are moved to the left of all the t a -matrices for the same nucleon, similar to (3): this way all the V 's will vanish when we square the amplitude.
Using the PV prescription to regulate the light-cone poles, the A-diagram contributions to the amplitude are The B-diagram contributions are The C diagrams are obtained by swapping the two projectile Wilson lines in Fig. 7, such that As explained in Appendix A, if one includes the contributions arising from moving the gluon exchanged between the two projectile quarks in the B-graphs across the cut, only the diagrams B 2 , B 9 and B 11 need to be calculated with the commutators of fundamental generators on the quark-1 line. Since B 2 = 0 we only need to take B 9 and B 11 from Eqs. (29i) and (29k) above and replace where the prime over the sum indicates that it is not a simple sum of the contributions in Eqs. (29), but a sum of all of those graphs and the contributions with the exchanged gluon moved across the cut. An expression similar to (31) can be obtained using Eq. (30) for the sum of C graphs (which would also include contributions arising from moving the gluons across the cut).
Adding all the A, B and C graphs by using Eqs. (28), (31) and (30) yields This is our final answer for the sum of A, B and C graphs in transverse momentum space.
Fourier transforming Eq. (32) into transverse coordinate space is straightforward. We use along with and Here, again, Λ ∼ Λ QCD is the IR cutoff in each individual nucleon.
Performing the Fourier transform of Eq. (32) we obtain This is our final answer for the A, B and C graphs. Previously a calculation of a similar quantity, the order-g 3 classical gluon production amplitude, was carried out in [52]. It would be desirable to compare the two results. Unfortunately a direct comparison appears to be impossible at the moment: first of all, the calculation in [52] appears to have employed a different regularization of α-integrals from our use of the b − 2 − b − 1 difference. The result obtained in [52] does not depend on b − 2 − b − 1 . In addition, the expression derived in [52] in momentum space still contains an analogue of our | k ⊥ × l ⊥ | in the denominator (see e.g. our Eq. (25) above), and is not simplified to cancel out this factor. It appears that significant simplifications of the result in [52] need to be performed before it could be compared with our Eq. (32). Furthermore, the calculation of [52] appears to be done in a gauge where there are no B or C graphs. Therefore, it is likely that the result of [52] does not involve the calculational tricks of moving the whole gluon across the cut in diagrams B and C detailed above around Fig. 4 and below in Appendix A. This could be another difference between the two calculations. Clearly more work is needed to compare our results with those in [52].

B. Graphs D and E
There is a subset of type-(i) diagrams (in the classification of Fig. 2) that we have not calculated yet. These are the diagrams where one of the nucleons in the projectile is a spectator in the amplitude, that is, it does not emit any schannel gluons but may interact with the shock wave, while another nucleon is a "spectator" in the complex-conjugate amplitude. An example of a diagram of this type representing the amplitude squared is shown in Fig. 9.
Clearly diagrams of the type shown to the left and right of the cut in Fig. 9 contribute at order-g 3 to the gluon production amplitude and need to be resummed. (Since both projectile valence quarks participate in the interaction in the amplitude squared, the contribution in Fig. 9 comes in with two powers of A C graphs.) We will label such diagrams as D-graphs: they are illustrated in Fig. 10. Let us stress that Fig. 10 only contains the D-diagrams which are not zero due to color averaging: we have to take a color trace (divided by N c ) for each projectile quark after squaring the amplitude. Taking into account that there are no emissions from quark 2 on the other side of the cut implies that many of the diagrams in this class give zero contribution to the amplitude squared: those graphs are not shown.
There are also E-diagrams, which are obtained by swapping the two projectile quarks (Wilson lines) in Fig. 10, Just like for the A, B and C graphs, in summing D-graphs one has to include the contribution where the gluons emitted by quark 2 are either emitter, absorbed, or both emitted and absorbed on either side of the cut. This again results in all gluon propagators being retarded. In addition, as explained in Appendix A, the sum of diagrams D 5 and D 6 , along with the diagrams D ′ 5 , D ′′ 5 , D ′ 6 , D ′′ 6 where either part of or the whole gluon loop are in the complex conjugate amplitude (see Fig. 15) is equal to diagram D 6 with the commutator of the color matrices [t a , t b ] instead of t a t b . With this in mind, the values of the D graphs are . (38e) Note also that color-averaging is implied: for instance, the part of the diagram D 1 contribution shown in Eq. (38a) is the one that survives a color-trace in the space of nucleon 2 after the amplitude is squared. More specifically when writing (t e t d ) 2 in Eqs. (38) we only keep the part of the expression which survives the (t e t d ) 2 → δ ed /(2N c ) substitution.
The sum of all the D-graphs (including the contributions from moving the gluons across the cut, as identified by the prime over the sum sign) is This result, along with the sum of all the E graphs, may also be directly obtained from Eq. (32) by "moving" the quark-gluon vertices on the line of quark 1 to the line of quark 2 in A, B and C graphs with an appropriate modification of color factors.
Fourier-transforming Eq. (39) we obtain This is our final result for the sum of the D graphs. Using Eq. (37) one can easily obtain the sum of the E graphs as well from Eq. (40). Again this result can also be obtained directly from Eq. (36).

C. Cross-checks
To test our main results for the order-g 3 gluon production amplitude in Eqs. (36) and (40) let us run a couple of cross checks.
First of all, if there is no shock wave, there should be no gluon production. A simple calculation shows that if one puts all U = 1 and all V = 1 then as expected. Another issue is gauge invariance. To test our results for gauge invariance, the order-g 3 gluon production amplitude can be calculated in different sub-gauges of the light-cone gauge, other than the PV sub-gauge in which the results in Eqs. (36) and (40) were obtained. In Appendix B we show that, while the calculations of the order-g 3 gluon production amplitude in the A ⊥ (x − → +∞) = 0 or A ⊥ (x − → −∞) = 0 sub-gauges are more involved than the PV sub-gauge calculation, the end result is the same in all three gauges for the unprimed sum of all A, B and C graphs, and for the unprimed sum of D and E graphs. Gauge-invariance also holds if we use retarded gluon "propagators". Note that obtaining commutators of fundamental SU(N c ) generators for the B and C graphs by moving gluons across the cut was shown to be a legitimate trick only in the PV sub-gauge of the light-cone gauge: it is likely that to re-derive the result from Eq. (36) in the A ⊥ (x − → +∞) = 0 and A ⊥ (x − → −∞) = 0 sub-gauges of the light-cone gauge, the gluon should also be moved across the cut in the diagrams like the first two graphs in Fig. 18, complicating the analysis. Hence, strictly-speaking, we have not shown that our result in Eq. (36) (with primed sums) is gauge-invariant: while it may still be gauge-invariant, at the moment we can think of it as a gauge-invariant amplitude (resulting from the unprimed sums) with the terms added (in the PV sub-gauge) casting it in a simpler form for obtaining the gluon production cross section.

IV. OUTLOOK
This paper is the first step in our project to calculate the first saturation correction in the projectile wave function to the p + A classical gluon production cross section. In this paper we calculated the order-g 3 gluon production amplitude. The main results are given in transverse coordinate space in Eqs. (36) and (40). Transverse momentum space amplitude is given in Eqs. (32) and (39). These results constitute the first saturation correction to the leadingorder gluon production amplitude in Eq. (3) and in Eq. (6) respectively.
To complete this project and find the order-α 3 s contribution to the classical gluon production cross section one needs to calculate the order-g 5 amplitude, again involving only two projectile nucleons (see Fig. 2). Similar to the D-graphs above, in calculating the order-g 5 amplitude we will be able to simplify the color algebra by employing the fact that only one of the two nucleons involved emits a gluon in the complex conjugate amplitude, as shown in the diagram (ii) of Fig. 2. Still, this appears to be a rather challenging calculation, and it is presently left for future work.

ACKNOWLEDGMENTS
The authors are grateful to Ian Balitsky for encouraging discussions. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics under Award Number DE-SC0004286.

Appendix A: Moving gluons across the cut
The purpose of this Appendix is to show that (a) the inclusion of the type-(iii) diagrams from Fig. 2 into the calculation can be easily accomplished by using retarded gluon Green functions instead of Feynman propagators; and that (b) the contributions of B and C graphs (shown in Fig. 7) to the O(g 3 ) amplitude along with the contributions of similar diagrams to the O(g 5 ) amplitude can all be efficiently found by calculating a subset of those graphs with color-commutators instead of the usual fundamental color factors. While these statements appear unrelated, they seem to be hard to disentangle for the B and C types of graphs.
Let us begin with contributions to the cross section coming from the squares of the A-graphs (defined in Fig. 6). It is straightforward to see by an explicit diagram analysis that the argument depicted in Fig. 3 applies in general, for the amplitude and the complex conjugate amplitude contributions coming from any two A-graphs. We illustrate this point in Fig. 11 by considering a slightly more involved example than that shown in Fig. 3, namely by studying the square of diagram A 1 . The argument of Eq. (11) applies here as well, after the following two observations. First of all, the interactions of the gluon that we moved across the cut in the second and third diagrams (to the left of the equal sign) of Fig. 11 with the shock wave cancel, as the gluon has the same transverse coordinate on both sides of the cut such that U U † = 1: hence we can treat this gluon as a free gluon propagator in the first graph, and as a free cut propagator in the second and third graphs. Second of all, in all three diagrams to the left of the equal sign of Fig. 11 the interactions of quark 1 with the shock wave cancel (for a similar reason to the gluon's, V V † = 1), leading to the same color factor on this quark line, tr [t a t b ] = δ ab /2, for all three graphs. Since the longitudinal coordinate integral ranges of the second and third graphs in Fig. 11, x + 1 > x + 2 and x + 2 > x + 1 , complement each other giving independent x + 1 and x + 2 integrals from −∞ to 0 each, we see that the sum of the three diagrams in Fig. 11 leads to a retarded gluon Green function denoted by an arrow in the last diagram in Fig. 11, just like it was demonstrated in Eq. (11). Other contributions to the amplitude squared made out of A graphs (both in the amplitude and in the complex conjugate amplitude) can be analyzed in a similar manner on a case-by-case basis, leading to the same result: all internal gluon lines should be represented by retarded Green functions. Note that to prove the use of retarded propagators for the gluon lines in the diagrams like A 2 one has to also employ the fact that the "+" component of the gluon momentum is conserved in the interaction with the shock wave. Now we have to analyze the contributions to the amplitude squared coming from B and C graphs, along with the interference terms between the A and B or C diagrams. The main tool here, in addition to Eq. (11), would be the cancellation of retarded Green functions from Fig. 5 (for the gluon exchanged between the quarks either before or after the shock wave both in the amplitude and the complex conjugate amplitude). Indeed one can proceed with the diagrammatic arguments, similar to how we outlined the analysis for the square of the sum of A diagrams above. While we performed this diagram-by-diagram check, it is rather lengthy and we are not going to present it here due to a very high number of diagrams one has to consider. Instead, a somewhat more compact approach to the problem is based on a formal argument which we present below.
Following [52] we denote gluon fields in the amplitude by A a µ , while the gluon fields in the complex conjugate amplitude are labeled by A a µ . These fields are not the background classical fields of the target, but rather "quantum" fields whose Wick contractions give us the s-channel gluon propagators used in calculating diagrams A, B, C, etc. By analogy to (5) we can define the Wilson line representing quark propagation over a finite x + -interval in the amplitude by The corresponding Wilson line in the complex conjugate amplitude would be V where A in the subscript denotes the fact that now one uses the field A a µ in the exponent of Eq. (A1).

Further we define gluon propagators as contractions
where D µν (k) is given by Eq. (16) while T and T denote time-and anti-time-ordering respectively. (As in the calculation leading to Eqs. (32) and (39) we will use the PV sub-gauge of the light-cone gauge.) Since only the upper "−" components of the gluon fields A a µ and A a µ enter the Wilson lines such as (A1), we will suppress the Lorentz indices below. That is, we will use Note that (cf. Eq. (11)) with D ab ret (x − y) denoting the retarded gluon Green function. Indeed we are setting up the well-known Keldysh-Schwinger formalism [57,58], which was used in small-x physics in [52,59,60].
Eqs. (A4) is a formal expression behind the cancellations like those shown in Fig. 5. To see this we concentrate on the Wilson lines representing quark propagators before the shock wave, both in the amplitude and in the complex conjugate amplitude. (For the cancellation it does not matter which interactions take place after the scattering in the shock wave.) The contribution of the diagrams in Fig. 5 is proportional to where the ⊗ sign indicates that the color indices of each product V † V are not contracted with each other, and are fixed at some arbitrary values. Time-ordering for A a µ -fields and anti-time-ordering for A a µ -fields are denoted by T A and T A correspondingly. For the classical gluon field diagrams we need contractions between the gluon fields A a µ and A a µ connecting the two projectile quarks: that is, we need contractions of the fields at ( b 1⊥ , b − 1 ) and the fields at . Indeed the diagrams in Fig. 5 are examples of such contractions. Expanding Eq. (A5) to the first order in the fields in each V † V , and analyzing only the b 1 b 2 contractions, we get This is the formal proof of the cancellation shown in Fig. 5. Note that before canceling, all the correlators assembled into the retarded Green functions. We have also employed an abbreviated notation by suppressing in the arguments of the correlators. A somewhat more involved calculation demonstrates that expanding Eq. (A5) to the second order in the fields in each V † V and concentrating on the b 1 b 2 contractions again, we would also get zero. This means that adding an extra gluon exchanged between the projectile quarks on either side of either diagram of Fig. 5 (as long as the gluon is both emitted and absorbed before the scattering in the shock wave) would still generate canceling diagrams. Two-gluon exchange is the highest order necessary in the classical field analysis [17]. Let us apply the Keldysh-Schwinger formalism described above to the analysis of the squares of B diagrams. We will work out the square of the diagrams B 9 + B 10 (as defined in Fig. 7), along with all the other similar graphs obtained from the square of B 9 + B 10 by moving gluons across the cut: contributions of other diagrams can be found by analogy. A sample of diagrams we need to sum up is given in the first line of Fig. 12. Namely, we fix the produced gluon to be emitted at light-cone time x + 1 in the amplitude and at x + 2 in the complex conjugate amplitude, and sum up all the graphs where the other two gluons are exchanged between the two projectile quarks, as long as the gluons are emitted and/or absorbed after the interaction with the shock wave in the amplitude and in the complex conjugate amplitude (that is, for x + > 0). Clearly the sample diagrams in the first line of Fig. 12 represent only a small subset of all the graphs that need to be included.
Instead of summing the diagrams we notice that due to the lack of gluon emissions and/or absorptions at x + < 0, the interactions of the projectile quarks with the shock wave cancel. Hence the contribution of all the B 9 + B 10 squared type of diagrams can be written as proportional to with the expression (A7) containing all the gluon exchanges between the projectile quarks. First let us rewrite Eq. (A7) as where the parenthesis around the V † V factor inside the first trace are placed there just to emphasize this term. Once again we are interested in contractions in Eq. (A8) connecting the b 1 and b 2 lines. By analogy to the analysis of Eq. (A5) one can conclude that two contractions between the fields in the second trace and the fields in the parenthesis inside the first trace cancel. We are left with contractions between , 0] A from the first trace and the same second trace, tr V † b2 [+∞, 0] A V b2 [+∞, 0] A . We may also have one contraction between the expression in the parenthesis of Eq. (A8) and the second trace, combined with another contraction either between , 0] A from the first trace and the second trace or between V † b1 [x + 1 , 0] A t a V b1 [x + 1 , 0] A from the first trace and the same second trace. This is exactly the answer illustrated in the second and third lines of Fig. 12. First of all, since the b 2 quark line brings in a factor of tr V † b2 [+∞, 0] A V b2 [+∞, 0] A , for each contraction of A − (b 1 ) coming from anywhere in the first trace in Eq. (A8) with A − (b 2 ) from the second trace there exists a contraction of A − (b 1 ) with A − (b 2 ); moreover, the second contraction comes in with a negative relative sign (due to hermitean conjugation in V † A as compared to V A ), completing the original contraction to a retarded Green function, in accordance with Eq. (A4a). Similarly, for each contraction between A − (b 1 ) with A − (b 2 ) there exists a contraction of A − (b 1 ) with A − (b 2 ), which comes in with the relative negative sign, again generating a retarded Green function, but now via Eq. (A4b). Let us stress again that since we are constructing gluon production cross section corresponding to the classical field, in this case we only need two contractions between the b 1 and b 2 lines, corresponding to the two gluons exchanged in the diagrams of Fig. 12. Hence our conclusions should be understood as valid for up to two gluon exchanges, but not necessarily beyond. The second point we need to make is that contractions with, say, , 0] A imply that the gluon in the amplitude can interact with the quark b 1 only at light-cone times 0 < x + < x + 1 , that is, to the left of the emitted gluon in the diagrams of Fig. 12. Moreover, since the interaction is identical to that with a gluon Wilson line connecting the points (0, b − 1 , b 1⊥ ) and (x + 1 , b − 1 , b 1⊥ ) along the x + light-cone direction. In passing we note that this is similar to the modification of the quark source current by a gluon exchange in the classical perturbative solution of Yang-Mills equations constructed in [12] (see Eqs. (8) there). To conclude we see that contractions between to the diagram to the left of the cut in the first panel of the second row of Fig. 12 (for one contraction) and to the second panel of the second row in Fig. 12 (for two contractions). The gluon propagators become retarded Green functions, and the color factors are replaced by commutators. Similar arguments for contractions between [+∞, 0] A give the right-hand side of the first diagram in the second line of Fig. 12 along with the last diagram in that line. (Let us stress that in the last two graphs in Fig. 12 we implicitly include a "crossed" contribution, where the quark-gluon vertices on quark line 2 are interchanged.) give the first diagram in the last row of Fig. 12 (more precisely, the part of the diagram to the left of the cut), with the last graph in this last row being a mirror image with respect to the final-state cut. Note that these last two diagrams in Fig. 12 are zero, since the contribution of quark-2 line is proportional to a color trace of a commutator. We see that to calculate the contribution to the cross section of the square of B 9 + B 10 + . . ., where the ellipsis represent the multitude of other graph in this class, we only need to calculate the square of B 9 with the retarded gluon Green function and with the color commutator instead of the standard fundamental color factor, along with the interference of the O(g 5 ) diagrams (containing a double commutator or two commutators, as shown in Fig. 12) with the leading-order (O(g)) gluon production amplitude, shown in the last four panels of Fig. 12. We have thus completed a demonstration of the original claim that to facilitate the calculation one can use retarded gluon propagators, calculate only a sub-set of B (and C) graphs with the commutators, along with calculating only the diagrams of the (i) and (ii) types by the classification of Fig. 2.
A complete proof of these results involves analyzing all the possible contributions to the amplitude squared. While still tedious, our formalism developed above makes the proof much more straightforward by reducing the vast number of diagrams one needs to consider to a much smaller number of Wilson line correlators like that in Eq. (A7). Let us illustrate our technique by a couple more examples. First consider the interference of the B 9 + B 10 + . . . type of graphs with C 5 (the upside-down reflection of B 5 from Fig. 7). This set of graphs is depicted in Fig. 13, where we only consider gluon exchanges between the projectile quarks at x + > 0 on the either side of the cut. Noticing that the shock wave interactions with the quark 1 cancel we conclude that the contribution of all the graphs of the type shown in Fig. 13 is proportional to We rewrite Eq. (A10) as If the two contractions are between V † b1 [+∞, 0] A V b1 [+∞, 0] A and V † b2 [+∞, 0] A V b2 [+∞, 0] A , then they cancel just like they did in the analysis of Eq. (A5). We are left with two options: one may have two contractions between , giving the answer shown in the first diagram after the equal sign in Fig. 13 Fig. 13. Interestingly the first diagram after the equal sign in Fig. 13 is zero, since the trace of the commutator we get in the color space of quark 1 is zero. Note that now the final answer in Fig. 13 is completely absorbed into the interference between the O(g 5 ) diagrams with the O(g) diagram: hence the diagram C 5 is absorbed into the O(g 5 ) gluon production amplitude, and is not included into our O(g 3 ) result presented in this work. Our final example demonstrating our formal technique involves the interference between the same B 9 + B 10 + . . . type of graphs with A 5 . It is illustrated in Fig. 14. Similar to the above we can see that the contribution of the diagrams from Fig. 14 is proportional to We rewrite this as An important difference now is that we are looking at only one contraction, corresponding to the single gluon exchange between the projectile quarks in Fig. 14. Similar to the above contractions between V † V 's in parenthesis cancel.
, 0] A in the second trace are zero: they are proportional to a trace of a commutator. We are left with contractions between is exactly the answer given in the last two diagrams of Fig. 14.
Note that one can repeat the above argument after moving the gluon line with color b (or c) in the first graph of The rest of the proof would be a more or less straightforward repetition of the above examples, which we verified but are not going to present here. In the end the propagators of all the internal gluon lines become retarded Green functions. The B and C graphs contribute as a subset of these graphs with color commutators.
This Appendix presents details of the correspondence between the "tree-level" Feynman diagrams and classical gluon fields. It is possible that an even more compact proof of our main results can be constructed: however, at the moment it is not known to the authors.  Fig. 10 along with the other graphs which need to be included. The emissions coming from quark 1 (see e.g. Fig. 9) in the complex conjugate amplitude are not shown as they are not relevant for the calculation at hand.
Our techniques presented above can be applied to the D and E graphs as well. However, in those cases a simple diagrammatic summation is possible. The D diagrams which need to be added together are shown in Fig. 15. (The emissions off of quark 1, such as those in Fig. 9, are not shown explicitly in Fig. 15 since they are not needed for our analysis.) In the end one obtains again in agreement with the main claims of this Appendix. The proof for the E graphs is obtained by a simple swap of the projectile quarks 1 and 2.

Appendix B: On Sub-Gauge Invariance
In this Appendix we demonstrate that up to order g 3 the gluon production amplitude does not depend on the particular sub-gauge used. Using we can write the gluon propagator in the A ⊥ (x − → +∞) = 0 sub-gauge of the light-cone gauge as Here we have split up the gluon propagator in the A ⊥ (x − → +∞) = 0 sub-gauge into two parts: one is equivalent to the PV gauge propagator, the other is proportional to δ(l + ). The gluon propagator in the A ⊥ (x − → −∞) = 0 sub-gauge can be written in a similar form, the only difference being the sign of the δ(l + ) term. The δ(l + ) component gives rise to extra terms for a given diagram when compared with the PV gauge. Due to gauge invariance of the amplitude it is expected that all these extra terms from all of the diagrams cancel out. This does end up being the case but the cancellations are not trivial. One needs to consider more diagrams than simply those shown in the main text as classes A, B, C, D and E.
There are three types of new contributions in the A ⊥ (x − → +∞) = 0 sub-gauge calculation of gluon production (as compared to the PV sub-gauge calculation): most of the A-graphs change (by ∆A i ) in the new gauge along with some of the D and E graphs (by ∆D i and ∆E i respectively), some of the B and C graphs acquire non-eikonal contributions which are still leading-order in the A ⊥ (x − → +∞) = 0 sub-gauge (we will refer to those contributions as "pinched") and there are new diagrams which can be identified as the shock-wave interaction corrections. All of these contributions will be described in detail below. These contributions are necessary when dealing with both the A ⊥ (x − → +∞) = 0 and A ⊥ (x − → −∞) = 0 sub-gauges of the light-cone gauge. Through the rest of this section we will be working in the A ⊥ (x − → +∞) = 0 sub-gauge. As will be shown, all of these new contributions can be written in terms of some sort of a "gauge rotation" acting on the order-g single gluon emission diagrams from Fig. 1. Since gauge invariance is valid separately for gluon emissions from quark 1 and quark 2, we will only perform the calculation explicitly for gluon emission from quark 1.
Pinched contributions originate from terms in the B and C diagrams which, while zero in the PV sub-gauge, are not zero in the A ⊥ (x − → +∞) = 0 sub-gauge. In order to see these extra contributions one cannot just treat the quark lines as Wilson lines, as was done in the rest of the paper; instead one must treat them as high energy quarks with large P + 1 and P + 2 momenta. Here we examine diagram B 3 depicted in detail in Fig. 16, to show how this pinched contribution arises. For the following calculation we only consider the left side of Fig. 16. That is, for the quark lines we only include terms up to the γ + from the shock wave. Here we show that this diagram can be written in terms of two components, the PV-sub-gauge contribution which comes from treating the quark lines as Wilson lines, and a "gauge rotation". The left (before the shock wave) side of the diagram gives Here b µ 1 = (0, b − 1 , b 1⊥ ) and b µ 2 = (0, b − 2 , b 2⊥ ) denote the intersection points of the projectile quarks trajectories with the shock wave. We also assume that the quarks are massless. We rewrite (B2) more compactly as where we have defined the following symbols: Before continuing with the remainder of the calculation it is useful to evaluate the following in the high-energy limit, Now let us focus on the shock-wave corrections. The diagrams which we have in mind are shown in the two left panels of Fig. 18 (cf. [41]). These diagrams do not belong to the A, B, C, D and E diagram types considered in the main text: the two-gluon interaction with the target nucleon in those graphs does not leave the nucleon in a color-singlet state, which is indicated by the nucleon breaking up in Fig. 18. One concludes that these diagrams are higher-order corrections to the Glauber-Mueller scattering in the target, which is limited by two gluons per nucleon [51]. However, his does not allow one to simply neglect these graphs, since this correction is enhanced by a power of A 1/3 1 (in the cross section) due to the extra projectile nucleon involved in the interaction. In the right panel of Fig. 18 we show a diagram in the A 1 class, with the interaction with the target limited to a single-gluon exchange with one of the target nucleons. Clearly all three diagrams in Fig. 18 are of the same order in our power counting: they are all order-g 5 and have the same projectile and target nucleons participating in the interaction. Hence the first two graphs in Fig. 18, along with other similar diagrams involving more multiple rescatterings in the target, are of the same order as the A, B, C, D and E diagrams and have to be included in the analysis. In fact one may worry why we did not include them in the main-text calculation in the PV sub-gauge of the light-cone gauge: below we justify neglecting these diagrams in the PV sub-gauge calculation performed in the main text.
To resum the corrections of the type shown in the left two diagrams of Fig. 18, we have to include corrections like this for either one of the many nucleons involved in the shock-wave interaction. We first consider scattering of two projectile quarks (coming from two different nucleons) on a shock-wave target. The corrections to the scattering on a single target nucleon are shown in Fig. 19 and labeled S 1 , S 2 and S 3 . (Our analysis will be similarly valid for scattering of the gluons coming from the projectile nucleons on the target.) Note that the diagram S 3 in Fig. 19 also has the standard eikonal contribution where the gluon exchanged between the projectile quark lines is long-lived in the s-channel: those types of contributions are included in the analysis of diagrams A − E in the main text and will be discarded here.
In the S diagrams shown in Fig. 19, quarks 1 and 2 are eikonal quarks traveling in the x + direction. 19. The diagrams correcting the interaction of two projectile quarks (top two straight horizontal lines) with a quark coming from a single nucleon in the shock wave (the thick red line).
labeled by a thick red line, originates in a nucleon from the target and thus has a large P − momentum and is traveling in the x − direction. While the evaluation of these diagrams is explicitly done for quarks, similar results exist for the case where quark(s) 1 and/or 2 are/is replaced by (a) gluon(s). For this calculation we will just focus on the quark case. The diagrams in Fig. 19 should be understood as being inside the shock wave; they can take place at any point in the shock wave. Since we are summing over all possible diagrams, in a given target charge distribution there will be diagrams where a given target nucleon has diagrams S 1 , S 2 and S 3 associated with it. We have to calculate and sum up all three diagrams associated with each nucleon in the target. As an example we calculate diagram S 1 explicitly.
FIG. 20. Diagram S1, a correction to the interaction with a single nucleon in the target shock-wave. The thick red line denotes a quark in the target nucleon.
In arriving at Eq. (B7) we have used the on-shell condition for the outgoing target quark, to eliminate the q + integral. Since quarks 1 and 2 are completely eikonal we replaced their spinor matrix elements by g α ′ − g β ′ − . In the target quark spinor matrix element from Eq. (B7) one can only have α = +, ⊥ and β = +, ⊥ since the gluon propagators vanish for either α = − or β = − due to the gauge choice. The leading contribution in P − 3 is given by the α, β =⊥ component. Naively one would expect such contribution to be sub-eikonal, suppressed by a power of P − 3 . As one could see below, this is indeed the case, but only until one integrates over l + picking the pole either at l + = 0 or at q + = 0. The residues at such poles generate enhancement (due to pinching of the poles in the product of propagators in (B7)), making the end contribution leading order. Since we are interested only in this pinched pole contribution, below we will evaluate all the expressions keeping the l + ≈ 0 approximation in mind.
Keeping the α, β =⊥ component allows us to write the product of the two gluon propagators in (B7) as Multiplying the spinor matrix element of quark 3 from (B7) by q ⊥α l ⊥β from the above expression yields The diagrams representing these four expressions are shown in Fig. 21. Notice that these satisfy the condition Let us go over the physical significance behind the four contributions in Eqs. (B14). The δ(x + 1 −x + 3 )δ(x + 2 −x + 3 ) term means all of these interactions are instantaneous, they occur at a single x + position. (In the end of the calculation all the x + coordinates are integrated out leading to the GM exponentiation of the result.) All of the diagrams in Fig. 21 are various corrections to the eikonal scattering seen in the MV model or GM approximation. Diagrams S ′ 1 and S ′ 2 consist of the eikonal scattering of the two quarks off the same target nucleon, each of them scattering in a classical field of the target quark, since the quark line is put on mass shell between the rescatterings (as indicated by a "cut" through the line). Indeed, each term in the square brackets of Eqs. (B14a) and (B14b) represents a single t-channel gluon exchange between the projectile quarks at b 1⊥ and b 2⊥ and the target quark at b 3⊥ . The contributions in (B14a) and (B14b) come in with a factor of 1/2 shown explicitly in the two left diagrams of Fig. 21. The only difference between S ′ 1 and S ′ 2 is the color factor associated with the target quark. Note that the contributions S ′ 1 and S ′ 2 also occur in the PV sub-gauge of the light-cone gauge. One can see this by noticing that they would remain if one drops the delta-function parts of the gluon propagators in Eq. (B7) (and in a similar calculation for S 2 ). Therefore, these contributions are sub-gauge invariant. In addition, these terms end up canceling out when one considers the amplitude squared (since bringing either of the t-channel gluons in these diagrams across the cut gives an overall minus sign [41]), meaning that these diagrams do not effect the gluon production cross section. This is why they are not included in the calculation in the main text leading to Eqs. (36) and (40). The other two contributions, S ′ 3 and S ′ 4 , which arise in the A ⊥ (x − → +∞) = 0 (and, up to a sign, in A ⊥ (x − → −∞) = 0) sub-gauge of the light-cone gauge but not in the PV sub-gauge, correspond to one of the projectile quark lines scattering in a classical field of the target while being color-rotated by the field of another projectile quark (see e.g. [17] for the color-rotation terminology). The first square brackets in each of Eqs. (B14c) and (B14d) contain a single t-channel gluon exchange corresponding to the classical field of the target quark. The second pair of square brackets in Eqs. (B14c) and (B14d) contain the "gauge rotation" illustrated by the dashed line in Fig. 21 using the notation defined in Eq. (B6) and Fig. 17. The arrow at the end of the dashed line indicates a commutator of fundamental color matrices for the dashed line "gluon" and for the "true" gluon entering the quark-gluon vertex: the first term in the commutator contains the color matrix for the dashed line placed to the left of the color matrix of the "true" gluon. To figure out how corrections from the diagrams S ′ 3 and S ′ 4 contribute to the amplitude with multiple rescatterings we need to add up all possible gauge rotations for multiple interactions with the target nucleons. The case of the two projectile quarks scattering on two target nucleons is shown in Fig. 22. On the left of the equality in Fig. 22 we have shown two contributions arising due to the corrections like those shown in the right two diagrams of Fig. 21: one for the first nucleons and one for the second nucleon. Analyzing the sum on the left we see that the gauge rotations inside the shock wave cancel out and we only end up with rotations on the outside of the charge distribution. The final result on the right of Fig. 23 gives an order-g 2 correction to the shock waves with the dashed line without an arrow defined as in Eq. (B6) and Fig. 17. Repeating this argument for any number of target nucleons, by including both corrections from the right two diagrams of Fig. 21 for each nucleon, we arrive at the same conclusion: the net result of all such corrections is equal to a diagram with a dashed line to the right of the shock wave minus the diagram with the dashed line to the left of the shock wave. We see that in the sub-gauge of interest each shock wave interaction is accompanied by the order-g 2 corrections as shown in Fig. 23. The dashed line contributions are given by the expression in the parenthesis of Eq. (B6), since all the x + coordinates are integrated out in going from S ′ 3 and S ′ 4 in Eq. (B14) to the contribution to the scattering amplitude.
Using this result it is straightforward to deduce the contributions of gauge rotations to the gluon production amplitude. One simply has to take an order-g gluon production amplitude for the scattering of two projectile quarks on the target and add all possible dashed lines to it (connecting a pair of s-channel lines) immediately to the left and to the right of the shock wave (with the appropriate signs). The resulting corrections to the single gluon production amplitude are depicted in Fig. 24   (II) shock wave corrections; (III) corrections to diagrams of type A and E. Each of the labels ∆Ai, ∆Bi, ∆Si, and ∆E2 denote both the diagram above it and the sign in front of the diagram. In the PV gauge all of these contributions are zero. Notice how the sum of all of these is zero, demonstrating sub-gauge invariance.
The last contributions we have to consider are the gauge dependent parts of the A, D and E diagrams. Calculating these is straightforward although time-consuming, just repeating the calculation done in the main text of the paper but this time including the part of the gluon propagators proportional to δ(l + ) (see Eq. (B1)). It is interesting to note that for an "instantaneous" term such as this, it does not matter whether one uses Feynman propagators or retarded gluon Green functions: hence our demonstration of the sub-gauge invariance also applies if one uses retarded gluon Green functions, as is done in the main text.
With this in mind all of the gauge-dependent terms of diagrams A and E (corresponding to a gluon emission from quark 1) are shown in Fig. 24 and are labeled ∆A i and ∆E 2 . These terms were obtained by an explicit calculation. Note that ∆A i = A i ( A ⊥ (x − → +∞) = 0 sub-gauge) − A i (PV sub-gauge), with the same definition for ∆E 2 . Also shown in Fig. 24 are "pinched" contributions ∆B i . Note that in Fig. 24 the labels ∆A i , ∆B i , ∆S i , and ∆E 2 include both the diagram above each of them and the sign in front of the diagram. All of the gauge-dependent contributions to Feynman diagrams at the order g 3 shown in Fig. 24, that is, the "pinched" contributions, shock wave corrections, and gauge-dependent diagram corrections cancel out, preserving gauge invariance. More precisely we have ∆A ′ 1 + ∆B 2 + ∆S 1 = 0; ∆A ′′ 1 + ∆S 2 = 0; ∆A 4 + ∆S 4 = 0; ∆A 5 + ∆B 10 + ∆S 7 = 0; ∆E 2 + ∆S 5 = 0; ∆B 5 + ∆S 3 = 0; ∆B 3 + ∆S 6 = 0.
Hence our result (the sum of all the A, B and C graphs and the sum of all the D and E graphs) is independent of the choice of sub-gauge in the light-cone gauge.