First Saturation Correction in High Energy Proton-Nucleus Collisions: II. Single Inclusive Semi-Hard Gluon Production

Exploiting recently obtained analytic solutions of classical Yang-Mills equations for higher order perturbations in the field of the dilute object (proton), we derive the complete first saturation correction to the single inclusive semi-hard gluon production in high energy proton-nucleus collisions by applying the Lehmann-Symanzik-Zimmermann reduction formula. We thus finalize the program started by Balitsky (see Ref.~\cite{Balitsky:2004rr}) and independently by Chirilli, Kovchegov and Wertepny (see Ref.~\cite{Chirilli:2015tea}) albeit using a very different approach to carry out our calculations. We extracted the functional dependence of gluon spectrum on the color charge densities of the colliding objects; thus our results can be used to evaluate complete first saturation correction to the double/multiple inclusive gluon productions.


Introduction
In Ref. [3], in the dilute-dense regime of the Color Glass Condensate (CGC) effective theory, we analytically solved the classical Yang-Mills equations to next-to-leading and nextto-next-to-leading orders in the field of the dilute object (the proton in proton-nucleus collisions). In this and follow-up papers, we turn our attention to computing observables; in particular, this paper is dedicated to using these solutions to derive a complete analytical expression for the first saturation correction to the single inclusive gluon production. This problem was previously considered in Refs. [1,2] for the single inclusive gluon production and in Refs. [4,5] for the double inclusive gluon production. However, these studies provided only incomplete results. Specifically, only order-g 3 single gluon production amplitude was computed in Ref. [2], while the order-g 5 production amplitude has never been computed in previous publications. In this paper, we give the complete first saturation correction to the single gluon production for the first time.
In the CGC framework, the leading contribution to particle production is due to classical gluon field. The quantum corrections, albeit important for certain aspects [6][7][8], are considered to be negligible in the first approximation. We follow this paradigm and work explicitly in the classical regime. It was shown that the classical approximation resums the powers of the saturation momenta of the projectile and the target [9]. The saturation momentum Q s is proportional to α s µ, where µ 2 is the hadron color charge squared per unit transverse area [10]. Therefore, the single inclusive gluon productions in high energy nuclear collisions can be expressed as a function of two dimensionless parameters Q 2 s,P /k 2 ⊥ and Q 2 s,T /k 2 ⊥ , involving the saturation scales of the projectile and the target respectively. In the notations of Ref. [2], (1.1) Here the overall factor of α −1 s manifests that we work in the classical approximation. The function f describes the most general situation of nucleus-nucleus scattering and is not known analytically. The series expansion of the function f in one of the arguments is termed the dilute-dense expansion. It is customary to consider the projectile as a dilute object and perform the expansion in powers of Q 2 s,P /k 2 ⊥ : (1. 2) The coefficients f n are analytically computable and can be non-perturbatively evaluated to all order in Q 2 s,T /k 2 ⊥ , at least in principle. In practice, however, only the leading term of this expansion was computed, the so-called dilute-dense approximation, see Refs. [11,12]: This truncation is justified as long as one is interested in the momentum region k ⊥ Q s,P of the gluon spectrum. As k ⊥ approaches Q s,P , higher order corrections have to be taken into account. For k ⊥ ≤ Q s,P , an all-order resummation is needed (see e.g. Refs. [13,14], and numerical simulations of Refs. [15][16][17][18]). We are mainly interested in the production of semi-hard gluons with transverse momentum k ⊥ > Q s,P . Thus our goal is to compute the next-to-leading correction f 2 in the expansion (1.2); the partial result for f 2 was presented in Ref. [2] where the contributions with amplitudes of order O(g 5 ) were not included. In Ref. [2], a diagrammatic approach based on the light-cone perturbation theory was used. In this paper, to extract single inclusive gluon production, we instead use the solutions of classical Yang-Mills equations derived in Ref. [3]. To this end we apply properly modified Lehmann-Symanzik-Zimmermann formula for Milne metric.
Our approach, besides the computational transparency, also allows us to go beyond configuration-averaged quantities; that is, we can perform "even-by-event" studies by evaluating particle production for a fixed distribution of the valence charge densities in both the projectile ρ P and the target ρ T . We aim to provide the functional dN d 2 kdy [ρ P , ρ T ]. The functional form of our result enables to study fluctuations of multiplicity; moreover having the functional form is important for two additional reasons. To appreciate the first reason, consider the configuration-averaged dN d 2 kdy ; it is often obtained from dN d 2 kdy [ρ P , ρ T ] by using the Gaussian McLerran-Venugopalan model [19,20]. This may miss some interesting effects potentially relevant for the observables. Indeed, we will demonstrate that, in general, the functional dN d 2 kdy [ρ P , ρ T ] is not an even function of ρ P , e.g. it contains a product of three color charge densities ρ P (x)ρ P (y)ρ P (z). The terms with odd powers of ρ P would vanish if averaged with a Gaussian distribution. However, for a proton projectile, a Gaussian model might not be the most appropriate in certain kinematic regions and a more elaborated approach might include non-trivial color charge correlations [21][22][23]. For instance, the proton model of Ref. [22] includes a nonzero C-odd part of the correlator of three color charges. Our resulting functional dN d 2 kdy [ρ P , ρ T ] will explicitly keep all even and odd powers of ρ P to the appropriate order and thus can be easily extended to account for non-zero three color charge correlations. The second reason is that the functional form of dN d 2 kdy [ρ P , ρ T ] allows us to evaluate multiparticle production in a very straightforward way [24][25][26] -by finding an average of the product of the functionals, e.g. for two-particle inclusive we simply get and analogously for n-gluon production. Although, the formulae provided in this paper can be easily applied to compute multigluon production; we will restrict this paper to evaluating the full functional including first saturation correction. We defer configuration-averaged gluon productions and correlations to a future work.
We also want to comment that our calculation properly includes interactions between multiple gluons at the same rapidity, e.g. the processes when the gluon from the proton wave function merges with the gluon emitted from the proton after the scattering off the target. The formalism used in Ref. [27,28] (see also reference therein for earlier works) does not account for such interactions.
Before moving on to the technical calculations, we would like to using Feynman diagrams to illustrate the first saturation corrections to single gluon production in high energy proton-nucleuse collisions. Fig. 1 shows diagram for the leading order single gluon production. The gluon could be emitted either before or after the interaction with the nucleus. Since there is only one gluon emission vertex and it only involves one color charge density, the amplitude is proportional to gρ b (x). It is important to point out that at the leading order each valence color charge independently emits gluons and scatters on the nucleus. They do not interact among themselves. Fig.2 shows part of the diagrams contributing to the order-g 3 gluon production amplitude. This is not meant to be exhaustive and there are other diagrams due to different vertex orderings and different light-cone time flow directions, see Ref. [2] for the complete ρ b (x) − Nucleus Figure 1: Schematic diagrams showing the leading order single gluon production in high energy proton-nucleus collisions. The ρ b (x) is the valence color charge density of the proton. The colored solid rectangular represents the highly Lorentz contracted nucleus. Figure 2: Partial diagrams representing the order-g 3 single gluon production amplitude. It involves interactions of two color charge densities ρ b 1 (x 1 ) and ρ b 2 (x 2 ).
set of diagrams. As can be seen, these diagrams contain three gluon vertices and involve interactions of two color charge densities, thus proportional to g 3 ρ b 1 (x 1 )ρ b 2 (x 2 ). Order-g 3 amplitude is the first part contributing to the first saturation correction to single gluon production. Not only the valence color charges could interact among themselves either before or after the collisions with the nucleus, the emitted gluons have final state interactions, as shown in the fourth and fifth diagrams in Fig. 2. Note that diagrams in Fig. 2 can be drawn from two copies of Fig. 1 with the appropriate attachments of gluon lines before or after collisions with the nucleus. Fig. 3 are six of the many diagrams contributing to order-g 5 gluon production amplitude. These diagrams contain five gluon vertices and they involve interactions of three color charge densities, thus proportinal to Figure 3: Partial diagrams showing the order-g 5 single gluon production amplitude. It involves interactions of three color charge densities be drawn from three copies of the order-g diagrams in Fig. 1 with the appropriate attachments of gluon lines. Both order-g 3 and order-g 5 amplitudes contribute to the first saturation correction to the single inclusive gluon production. One can symbolically represent the single gluon production amplitude M = M (1) + M (3) + M (5) + . . .. Here the subscript indicates the power of coupling constant g. However, it should be pointed out that this is not exactly a perturbative expansion in terms of coupling constant because we only take into account diagrams that are enhanced by the color charge density of the proton at each order. Therefore M (1) ∼ gρ P , M (3) ∼ g 3 ρ 2 P and M (5) ∼ g 5 ρ 3 P . The single gluon production is proportional to the amplitude squared The crossing terms between order-g and order-g 3 amplitudes are proportional to odd powers of the color charge density In McLerran-Venugopalan model (see also Ref. [29]), ρ P satisfies a Gaussian distribution; therefore ρ 3 P terms vanish upon averaging over all possible configurations of color charge density. The first saturation correction therefore starts from Both the order-g 3 amplitude squared and the crossing terms between order-g and order-g 5 terms are proportional to g 6 ρ 4 P . In a Gaussian model for the color charge density, one has g 4 ρ 2 P ∼ g 4 µ 2 P ∼ Q 2 s,P so that the first saturation correction is proportional to 1 αs Q 4 s,P .
This should be compared with the leading order term |M (1) | 2 ∼ 1 αs Q 2 s,P . This power counting in terms of the saturation scale has already been formally shown in Eq. (1.2). The goal of the paper is to derive M (3) and M (5) and finish the computation of the first saturation correction to single gluon production (previous studies in the literature only computed M (3) ).
The paper is organized as follows. In Sec. 2, we briefly review the key results of the previous paper [3]. In particular, we set up an initial value problem for the classical field in the forward light cone and list the required solutions. In Sec. 3, in order to extract the single inclusive gluon production, we detail the Lehmann-Symanzik-Zimmermann reduction formula in the Milne coordinate system. We finally turn to genuinely novel material in Secs. 4 and 5, where for the first time we compute the functional dN d 2 kdy [ρ P , ρ T ] including the first saturation correction in the projectile field. We finish with conclusions in Sec. 6. For the sake of completeness, we provide more technical details in the appendices.

General Setup and Solutions of Classical Yang-Mills Equations
We briefly recapitulate the general setup of our calculations and the results obtained in Ref. [3]. In the CGC framework, the classical gluon fields produced in high energy proton-nucleus collisions are determined by the classical Yang-Mills equations D µ F µν = J ν . Here the color charge current for a right-moving projectile and a left-moving target In the forward light cone x + > 0 and x − > 0, the color charge current vanishes and the effects of the sources ρ P , ρ T are reflected in the initial conditions. We adopt the procedure described in Refs. [12,30] and consider the Fock-Schwinger gauge x − A + + x + A − = 0. This gauge suggests a convenient parameterization for the solutions, which are explicitly assumed to be boost invariant, In terms of α(τ, x), α i (τ, x), the classical Yang-Mills equations read with the initial conditions [30,31]: The α i P (x) and α i T (x) are the Weizsacker-Williams (WW) gluon fields of the projectile and target, respectively; they depend on the transverse coordinate and are related to ρ P (T ) through the static Yang-Mills equations ∂ i α i P (T ) (x) = gρ P (T ) . Since each field α i P (T ) (x) is a pure gauge, it is convenient to use the following parametrization in terms of Wilson lines V for the projectile and U for the target α i In general, the WW field of the projectile (target) contains all powers in ρ P (ρ T ); for our purpose only a power series expansion in ρ P is required: 3) The first three expansion coefficients are where the gluon exchange potential with m being some appropriate IR cutoff.
The initial value problem can be solved order-by-order in ρ P while preserving nonperturbative dependence to all orders in ρ T . In order to facilitate this, we consider a subgauge transformation defined by the condition ∂ i β i (τ = 0, x) = 0 to which we will refer to as the initial time Coulomb gauge, see e.g. Ref. [32,33]. The fields β and β i are non-linearly related to α and α i through a gauge transformation (2.8) The Ω = Ω(x) is fixed to satisfy ∂ i β i (τ = 0, x) = 0. A closed form of Ω is not known but a perturbative expansion in coupling constant is obtained in Ref. [3]. The subgauge transformations, not affecting the FS gauge, preserve the form of the Yang-Mills equations and only modify the initial conditions. Below we will present solutions of the classical Yang-Mills equations as power series expansion in the coupling constant g in the initial time Coulomb subgauge (keeping all power of gρ T ). Due to the structure of the classical Yang-Mills equations, only odd powers of g are non-vanishing, that is The order-g and order-g 3 solutions in momentum space are Here (2.14) The above solutions are expressed in terms of order-g initial conditions b η (k), b ⊥ (k) and order-g 3 initial conditions β (3) (τ = 0, q), β ⊥ (τ = 0, q). To present expressions for the initial conditions, we introduce an auxiliary quantity which will repeatedly appear in our calculations: The physical meaning of this combination is clear: it represents the eikonally color rotated projectile gluon field by the target Wilson line. We will also use the expansion of ζ i (x) in terms of coupling constant g with the coefficients For example, the explicit expression for ζ i (x) at order-g is By analogy with ζ i (x), we also introduce the eikonally rotated projectile color charge density In terms of these auxiliary quantities, we have order-g initial conditions b η (k) = 2β (1) (2. 19) In order to make the further notation uniform we also introduce In other words we can expand the two dimensional vector ζ (1),i into the longitudinal and transverse parts The order-g 3 initial conditions are Although we will not need the order-g 5 solutions of the classical Yang-Mills equations when computing order-g 5 gluon production amplitude using the LSZ reduction formula, we do need the order-g 5 initial conditions. The initial conditions at order-g 5 are β (5) (2.24) (2.25)

Defining the Semi-Hard Gluon Productions
In this section we conform the conventional four-dimensional Minkowski LSZ procedure (see the textbook [34]) to the problem at hand, that is 2+1 Milne coordinate system with the initial conditions defined at τ = 0. There are two main sources of apparent differences from the conventional LSZ: the first is due to a non-trivial metric modifying the free-field modes and the second is due to a non-vanishing field at the initial hypersurface τ = 0.

The LSZ reduction formula in Milne coordinate
We begin by constructing free-field modes for the dynamical degrees of freedom. In our case, for the gluon field there are two dynamical degrees of freedom; they are represented by two fields β(τ, x) and β ⊥ (τ, x). The remaining field Λ(τ, x) is non-dynamical and does not contribute to the particle spectrum. The free-field equations forβ(τ, x) ≡ τ β(τ, x) and β ⊥ (τ, x) are The free-field solutions can be easily found in the momentum space. We choose the solutions in a such a way as to resemble the plane-wave behaviour at asymptotically large τ . This fixes the form of the functions -they are Hankel functions of first and second kind (see Ref. [35] for mathematical definition); the first (zeroes) order Hankel function solves the first (second) equation in Eqs. (3.1). The normalization can be fixed by requiring the consistency between the equal time commutation relations for the fields and the canonical commutation relations for creation and annihilation operators, as we demonstrate below. We thus have two sets of mode functions and They satisfy the Wronskian identities Next we construct the mode expansions for both dynamical fields in terms of creation and annihilation operators (3.5) One can check that the normalization of the mode functions leads to the consistent commutation relations, that is for the fields and for the creation and annihilation operators. The mode expansion (3.5) can be inverted in order to express the creation operators in terms of the fieldŝ (3.8) Here . Similar to Minkovski spacetime, we can start from the trivial identity for the creation operator for the asymptotic out state (3.10) Here we used the fact that H (H (2) 0 ) satisfies the first (second) free field equation in Eqs. (3.1). Additionally we took into account that the fieldsβ(τ, p) and β ⊥ (τ, p) are the (formal) solutions of the non-linear Yang-Mills equations with the source terms S η (τ, p) and S ⊥ (τ, p). Here we see that the creation operators for the asymptotic out state have two distinct contributions. The first is due to the initial "surface" contribution τ = 0, while the second due to the "bulk" evolution in the forward light cone.
We can also computeâ † p (τ = 0) andĉ † p (τ = 0): (3.12) Where the constants C 1 and D 1 can be found using the initial conditions β(τ = 0, p) = p ⊥ 2 C 1 and β ⊥ (τ = 0, p) = D 1 . Following the notations of Ref. [4] we define the surface contributions S η and S ⊥ by (3.13) and the bulk terms B η and B ⊥ by (3.14) Thus the LSZ reduction formula can be written in the compact form As the reader will see below, differentiating between the surface (fully defined by the initial conditions at τ = 0) and the bulk terms (determined by the final state interactions) helps to pin-point the importance of the final state interactions.

Single inclusive gluon spectrum
Using the creation and annihilation we define the single inclusive gluon spectrum (for the fixed valence charge densities) as (3.16) Here γ denotes two gluon polarizations γ = η, ⊥; we sum over the polarizations as we are only interested in the polarization-blind single inclusive gluon production. It is instructive to analyze the k ↔ −k symmetry property of each term in dN/d 2 kdy as it does not require explicit expressions. This analysis will be helpful in isolating the odd harmonics in double gluon productions. As a first step, let's discuss the symmetries of the surface and bulk contributions.
The initial WW fields are real functions of x; consequently the surface terms have the following symmetry in the momentum space The bulk terms require a more detailed analysis. Superficially, one can argue that since the source terms S a η (τ, x), S a ⊥ (τ, x) are real functions of x and consequently S a, * η(⊥) (τ, k) = S a η(⊥) (τ, −k), we should expect that the entire bulk contribution has to be symmetric under k ↔ −k. This argument however misses the fact that the Hankel functions in the definition (3.14) are complex valued. Taking this into account and explicitly separating the Hankel functions into the real and imaginary parts we can write Note that in the configuration space, the first term is purely imaginary while the second term is purely real. There is an overall phase difference between two term of π/2, see a related discussion in Ref. [5]. The newly introduced functions satisfy Similar separation can be defined for the transverse component of the bulk term with the properties We are now in position to discuss the symmetry properties of each term in Eq. (3.16). The pure surface term is symmetric with respect to k ↔ −k.
The interference terms between the surface and bulk contributions have two contributions with different properties under k ↔ −k (3.23) We grouped the terms according to the symmetry: the first (second) piece involving B a γ (B a γ ) is antisymmetric (symmetric) with respect to k ↔ −k. As was identified in Ref. [4], the antisymmetric part is ultimately responsible for generating a non-vanishing contribution to odd harmonics at the lowest non-trivial order. The pure bulk term reads The contribution in the first bracket is symmetric while in the second bracket is antisymmetric under k ↔ −k. To summarize this discussion, the bulk terms are complex-values and their imaginary parts contribute to odd harmonics either through the bulk surface interference or the pure bulk terms. Note that overall dN d 2 kdy is real. In practice, our calculations will be done for power series expansion for the single gluon distribution. We have This power series expansions were terminate at the fifth order, as it is sufficient to extract the first saturation contribution to single inclusive gluon production, which is up to order g 6 : In the following we calculate each order n (i) (k) separately. From our previous analysis, it is obvious that, for a fixed configuration of the valence charges, the leading order n (2) (k) is even under k ↔ −k while all higher order terms have both even and odd contributions. Note that g 4 terms is cubic in ρ P and thus, at least in the Gaussian MV model, does not contribute to the configuration/event averaged single inclusive gluon production.

The Gluon Production Amplitude
We now proceed with assembling all the pieces together in order to find the functional dN d 2 kdy [ρ P , ρ T ]. First we explicitly write down the surface and the bulk terms.

The surface terms
The surface terms are defined by the initial condition on the light cone τ = 0; thus to find them one does not require solving the classical Yang-Mills equations in the forward light cone τ > 0. Starting from the definitions we obtain that the order-g surface terms are simply The order-g 3 and order-g 5 contribution to the surface terms can be also extracted straightforwardly (4.6)

The bulk terms
In this section we evaluate the bulk terms B η,⊥ and B The bulk terms B η and B (3) ⊥ The order-g 3 bulk terms are defined by Here the corresponding source terms are From Eq. (4.8) and Eq. (4.9) we obtain the order-g 3 bulk terms (4.10) In computing the bulk terms, we need to calculate the proper time integral from τ = 0 to τ = ∞. The integrands involve products of Bessel functions. The integral formulas used are summarized in Appendix A. Results in eq. (4.10) have already been obtained in ref. [4]. We now move on to compute the new order-g 5 bulk terms.
The bulk term B (5) η (k) The first task is to compute explicitly the source term at order-g 5 We substitute the solutions i (τ, k) given in Eqs. (2.10), (2.11), (2.12), and (2.13) in order to obtain the explicit expression for S (5) η (τ, k), (4.12) Using the above result, we can evaluate the bulk term B Again, it requires integrations over the proper time from τ = 0 to τ = ∞. All the integrals involved have closed form results. The integral formulas used are given in the Appendix A. We only present the final result here (4.14) Here I 1,2,3 (p, q, k) are functions of the momentum only and do not depend on the target Wilson line or the projectile color charge density. The introduction of these three functions is because the angular integrations cannot be carried out analytically. We therefore denote the various parts containing the angular integrals by these auxilliary functions: The definition of the function L 011 (a, b, c) can be found in the Appendix A. It is expressed in terms of simple elementary functions although containing Heaviside step functions depending on the relative magnitude of a, b, c. Superficially it appear that functions I 1,2,3 (p, q, k) might have multiple singularities due to vanishing denominators in the integrands. We demonstrate in Appendix C that all singularities are removable.
The bulk term B For the transverse polarization, the corresponding order-g 5 source term has the expression Again, using the solutions β (1) i (τ, k) given in Eqs. (2.10), (2.11) (2.12), and (2.13), one obtains its explicit expression The source term is then used to compute the order-g 5 bulk term, contributing to single gluon production amplitude, by We skip the details of doing the proper time integral and present the final result here (4.21) Similar to the analysis for B (5) η (k) in the previous section, we introduced another three auxilliary functions G 1,2,3 (p, q, k) defined as (4.24) The explicit expression of functions L 000 (a, b, c) and L 110 (a, b, c) are given in Appendix A. The analytic structure of these functions are studied in Appendix C. Althought it is impossible to obtain closed form expressions for these auxilliary functions, after taking care of the removable singularities, they are easily evaluated numerically for given momenta p, q, k.

Calculating the Single Inclusive Gluon Production
We are now ready to proceed with evaluating the single gluon production up to order g 6 . We will eventually express the gluon production in terms of the proton color charge density ρ a P (k) and the target Wilson line U ab (k). The following explicit expressions will be repeatedly used. We list them for reference.
Further more, the eikonally color rotated order-g 3 and order-g 5 Weiszacker-Williams fields are

Leading order gluon production
The leading order gluon production is defined by Using the explicit forms for we obtain that at the leading order This result reproduces the dilute-dense single gluon spectrum calculated in Refs. [12,36,37]. Note that this expression is given in the functional form with explicit dependence on ρ P and ρ T . For comparison with high order corrections, we use the Lipatov vertex [38][39][40] Then the leading order gluon production can be expressed in a more informative form The color structure can be shown by Fig. 4. The coefficient function is just the Lipatov vertex squared.

Next-to-leading order gluon production
The order-g 4 contribution to the single gluon production is defined by The surface terms S ⊥ (k) and S Collecting the results for the surface and bulk terms at order-g 3 , we have (5.13) From these surface terms and bulk terms, after some tediou algebra, substituting eqs (5.1), one obtains the next-to-leading order single inclusive gluon production Here c.c. represents the complex conjugate terms and we have introduced the shorthand notation p = d 2 p (2π) 2 . The color structure (functional dependence on ρ P and U ) are correctly captured by the two diagrams in Fig. 5. The kinematic function E 1 has the expression One immediately sees that the complex conjugate amplitude is represented by the Lipatov vertex We thus introduce an effective vertex at order-g 3 , (5.17) The coefficient function in eq. (5.15) follows simply as the product of the two effective vertices.
See Fig. 5a for the momentum assignments and the color structure.
The second kinematic function in eq. (5.14) is Again, isolating the Lipatov vertex in the complex conjugate amplitude we introduce another effective vertex at order-g 3 , As a consequence, the coefficient function E 2 can be expressed as a product of the two effective vertices Fig. 5b for the assignments of the various momenta.
Note that E 2 has both real and imaginary parts. The imaginary part is always accompanied by the sign function sgn(k × p) = (k × p)/|k × p|. It is precisely this imaginary part that will be responsible for the asymmetry between n (4) (k) and n (4) (−k). In n (4) (k), the odd part under k ↔ −k has already been extracted in Ref. [4]. The results were contrasted with those obtained from a diagrammatic approach of Ref. [2] in Ref. [5]. The even part under k ↔ −k is computed here for the first time. Also note that since n (4) (k) is cubic in ρ P , it does not result in any nontrivial contribution to the single inclusive gluon production after averaging over the projectile color charge densities in the McLerran-Venugopalan model. However, it will be important for double inclusive gluon productions (as n (4) (k)n (4) (p) is proportional to ρ 6 P , which does not vanish in the McLerran-Venugopalan model).
Considering the color structures for n (4) (k) in Eq. (5.14), the terms can be grouped into two categories according to the number of target Wilson lines involved: two or three Wilson lines. From Fig. 1 and Fig. 2, there seem to be a lot of possible combinations of diagrams when multiplying the order-g 3 amplitude and the order-g complex conjugate amplitude. However, as far as color and transverse momentum flows are concerned, all possible graphs can be properly represented by the two diagrams in Fig. 5. With the help of the effective vertices L j (p 3 , k), Γ j 1 (p, p 1 ), Γ j 2 (p, p 1 , p 2 , k) , one only needs to directly compute the two diagrams shown in fig. 5 to obtain the correct next to leading order correction to the single gluon production.
Finally, as a consistence check, when the target Wilson line U (x) = 1, meaning the absence of the nucleus, there will be no gluon productions. Using U ab (p) → δ ab δ (2) (p), it is easy to check that n (4) (k) expressed in Eq. (5.14) vanishes.

First saturation correction to single gluon production
In this section, we calculate the order-g 6 contribution to the single inclusive gluon production, the so-called first saturation correction. This contribution is quartic in ρ P and will give the first non-vanishing correction to the configuration-averaged single inclusive gluon production in the McLerran-Venugopalan model.
At order-g 6 , the single gluon production is defined by

(5.23)
We have obtained all the surface and bulk terms needed to compute the final expression for n (6) (k) (see Sects. 4.1 and 4.2). The explicit calculation is very long and tedious (it needs substitution of eqs (5.1)). We thus only report the final result here. All the terms are grouped into three categories according to the number of target Wilson lines they contain: two, three and four Wilson lines.
We start with the order-g 3 amplitude squared The three scalar functions H 1 (p, q, p 2 , p 4 ), H 2 (p, q, p 2 , p 3 , p 4 , k), H 3 (p, q, p 1 , p 2 , p 3 , p 4 , k) only depend on transverse momentum. We present their explicit expressions in the following. (5.25) We have used the identity (p 2 × q)p 2 ⊥ − (p × q)(p · p 2 ) = (p · q)(p 2 × p). It is obvious that H 1 is a real function of transverse momentum. Up to a minus sign, H 1 is nothing but the product of the two effective vertex introduced in Eq. (5.17) Fig. 6a for the momentum assignments.
There are six terms in H 2 ,
There are thirteen terms in H 3 , It is an easy exercise to show that H 3 can be expressed in terms of the effective vertex we introduced in Eq. (5.21) Fig. 6c for the various momentum assignments in the amplitude and the complex conjugate amplitude.
There are two sign functions present in the expression of H 3 , sgn(k×p) and sgn(k×q). Imginary parts of H 3 contain only one of the sign functions but not two at the same time. On the other hand, some terms in the real part contain product of both sign functions. It is clear that a single sign function is responsible for breaking the symmetry between k ↔ −k.
As for the color structure of the order-g 3 amplitude squared Eq. (5.24), although there are many diagrams when multiplying the amplitude with the conjugate amplitude, according to the number of Wilson lines, they can be representd by the three diagrams in Fig. 6. Depending on how many gluon lines crossing the target nucleus in the diagram, they corresponds to two-, three-and four-Wilson line terms. With the help of the effective vertices Eqs. (5.17) and (5.21), the precise expression can be computed directly from these three diagrams.
We move on to present the final result for the crossing terms between the order-g and order-g 5 gluon production amplitudes. Again, the terms are grouped into three categories according to the number of target Wilson lines.
The three kinematic functions F 1 , F 2 , F 3 only depend on transverse momentum. Their explicit expressions are presented in the following. For F 1 , (5.32) Subtracting the Lipatov vertex in the complex conjugate amplitude we introduce an effective vertex at order-g 5 , Then the coefficient function can be expressed as a product of two effective vertices Fig. 7a for the momentum assignments. It is interesting to note that the effective vertex Υ j 1 (p, q, q 2 ) is independent of the momentum k. It is also apparent that F 1 is a real function.
There are eight terms in F 2 . The final result can also be expressed as a product of two effective vertices We decompose the order-g 5 effective vertex Υ j 2 into two parts Υ j 2 (p, q, p 2 , p 3 , k) = Υ 2 (p, q, p 2 , p 3 , k) The two functions are See Fig. 7a for the momentum assignments. F 2 contains both real part and imaginary parts. The imaginary part is associated with the sign function sgn(k × q).
There are nineteen terms in F 3 . Six of the terms involve the auxiliary functions I 1,2,3 and G 1,2,3 we introduced in Sec. 4.2. These functions contain definite angular integrals that no closed form results are known to us. The computations are straightforward but very tedious. The final expression of F 3 can also be expressed as a product of two effective vertices The third effective vertex at order-g 5 , Υ j 3 , can be further decomposed into two parts depending whether it is parallel or perpendicular to the momentum k, Here the two functions Υ 3 and Υ ⊥ 3 are saturation correction to single gluon production There are a few comments about n (6) (k). First of all, it is impossible to relabel the integration momenta to make the color structures of the terms that have equal number of Wilson lines the same. Taking the two terms having four Wilson lines as an example, denoting the eikonally rotated color charge density ρ b R (k) = ρ a P (p)U ab (k−p), the two color structures are

(5.45)
Here k is some given external momentum. The integration variables are p, q. It is obvious that one can not redefine p, q to make these two expressions the same. On the other hand, if one is interested in the total number of gluons rather than the gluon spectrum, with an additional integration d 2 k in the game, one can redefine q → k and k → −q in the second expression so that the two color structures are the same. There is another possibility that after averaging over the projectile color charge configuration, terms with equal number of Wilson lines might be combined. We have explicitly checked this situation using the MV model, however, only partial terms could be combined and no significant cancellations happen.
The n (6) (k) is expressed as a functional of the color charge densities ρ P and ρ T (through the Wilson line U). This event-by-event expression is useful for computing double/multiple gluon productions. To obtain event averaged gluon spectrum one needs to perform ensemble average over all possible color charge configurations of the projectile and the target. One popular statistical distribution for the color charge density used in the literature is the Gaussian distribution. It is straightforward to carry out the average over ρ P . However, the average over ρ T contains four adjoint Wilson lines and is in general very complicated albeit possible [41][42][43]. There also exists other approximation schemes, see discussions in Refs. [44][45][46]. Event averaging will be studied in a follow-up paper.
As far as color structure is concerned, n (6) (k) can be simply computed from the six diagrams in Fig. 6 and Fig. 7. These six diagrams capture the correct dependence on ρ P and ρ T . Of course, the realistic coefficients F 1,2,3 and H 1,2,3 can not be read from the diagrams. However, with the help of the effective vertices, one can directly computing the correct expression for n (6) (k) from the six diagrams .
Compared with the leading order gluon production n (2) (k), the most important feature of the first saturation correction n (6) (k) is that it includes nontrivial final state interactions. This is explicitly shown in figs. 6b, 6c, 7b and 7c. Gluons produced by different color sources merge together after independently scattering on the target nucleus. Figs. 6a and 7a only contain initial state effects, no gluon interactions after scattering with the nucleus.
Finally, as a consistence check, when the target Wilson line U (x) = 1, Eq. (5.44) vanishes as expected because no gluon will be produced if there are no interactions with the nucleus.

Conclusions
In this tour de force calculations, we obtained the first saturation correction to the single inclusive semi-hard gluon production in high energy proton-nucleus collisions. The main result is given in Eq. (5.44). It is expressed as a functional of the color charge densities of the projectile ρ p and the target ρ T . Interestingly we found that to obtain the correct functinal dependence on ρ P (T ) , we only need to compute the six diagrams shown in Fig. 6 and Fig. 7. The coefficients of these functional dependence, however, have to be determined by a full analysis of all the contributing diagrams or through the tedious computation using the LSZ reduction formula, as have been done in this paper. To encapsulate our results, we introduced several effective vertices. With these effective vertices, directly computing the six diagrams in Fig. 6 and Fig. 7 leads to the correct result for the first saturation correction to single gluon production.
The order-g effective vertex L j is the famous Lipatov vertex. For the two effective vertices at order-g 3 , Γ j 1 , Γ j 2 and the three effective vertices at order-g 5 , Υ j 1 , Υ j 2 , Υ j 3 , although they are very convenient quantities to organize the calculations, they represent summations of many diagrams and it is hard to percieve their physical meanings.
Eq. (5.44) is an event-by-event result, in a follow-up paper, we will averaging over the color configurations. The averaging heavily depends on different modelings of the statistical distributions of ρ P and ρ T . It is probably more reasonable to model the dilute proton and the dense nucleus differently.
Our ultimate goal is to evaluate the saturation effects of the dilute object on double gluon productions in high energy dilute-dense scatterings, this includes proton-nucleus collisions as well as light-heavy ion collisions. Since the first saturation effects also contains non trivial final state interactions, we hope to study this effect, work out the phenomenology and compare with data from LHC and RHIC (see early works [47,48]).
In appendix. D, we made comparisons with the order-g 3 gluon production amplitude derived using the light-cone perturbation theory in [2]. The odd part under k ↔ −k of the order-g 3 amplitude is shown to be exactly the same even though their calculation was done in the light-cone gauge while ours is done in the Fock-Schwinger gauge. Other than that, comparisons at amplitude level at two different gauges would not make sense. Only the gauge invariant gluon production should be compared.
It is also possible that beyond the leading order, rapidity dependent fluctuations might play an important role in semi-hard gluon productions, see Refs. [6,7] for their influences on the bulk properties of the Glasma (recall that we assumed boost invariance in this paper). It would be interesting to evaluates its contribution compared with those from the first saturation correction. In this sense, our complete result on the first saturation correction, besides its own importance, also serves as a very useful reference.

A Definite integrals involving Bessel functions
In order to simplify the readers task to reproduce our results, we list all definite integrals involving Bessel function that were used in the main body of this paper.

Integrals involving one Bessel function
Here we used the following regularization procedure in order to assign the meaning to the integrals over J 0 and Y 0 : • Taking the limit p ⊥ → 0, in the closure relation (A.4) • Additionally from Eq. (A.6), we have

Integrals involving two Bessel functions
The closure relation for Bessel functions is Together with the crossproducts (A.6) the closure relation leads to

Integrals involving three Bessel functions
Here we list the most general expressions; they are followed by the specific integrals used in the main body of the paper.

General formulae
The following integrals contain products of three Bessel functions in the integrands. (see more in Ref. [49]): • Integrand xJ 0 (ax)J 1 (bx)H 1 (cx) Lets consider the real and imaginary part separately, H is non-zero only if |c − b| < a < c + b.
The imaginary part of the integrand involves Y 1 (cx). The associated integral is given by In order to simplify notation it is convenient to introduce (A.12) In particular, for |c − b| < a < c + b, the above expression simplifies to 0 (cx) For the real part, we quote the formula from [49] (A.14) It is nonvanishing only in the parameter region |a − b| < c < a + b. Here and a, b, c > 0.
It is handy to use the notation for the explicit form In particular, when |a − b| < c < a + b, it simplifies to (A.32)

Specific formulae
Using the above general forms for the integrals involving three Bessel functions, it is straightforward to establish the following list of identities

Integrals involving four Bessel Functions
We also need integrals involving four Bessel Functions. Our strategy is to reduce the number of Bessel functions in the integrand from four to three, which we have explicit formula do the integration.
0 (k ⊥ τ ): By applying the Garf's formula (see Appendix of Ref. [3]) and we perform the desired reduction to get With this reduction we obtain In this expression, one can substitute for e iΨ , since only even part of e iΨ contributes to the final result. The same applies to Eq. (A.47).
B Relations among L 000 (a, b, c), L 110 (a, b, c) and L 011 (a, b, c) We previously introduced three functions L 000 (a, b, c), L 110 (a, b, c) and L 011 (a, b, c). It may appear that they are independent. However, further analysis shows that only one function is independent. Indeed, using the reduction (see Appendix of Ref. [3]) cos φJ 0 (wτ ) (B.1) with w 2 = a 2 + b 2 − 2ab cos φ and Eq. (A.8), we get Note that the final expression is symmetric with respect to a ↔ b. This establishes the connection between L 110 (a, b, c) and L 000 (a, b, c).
For L 011 (a, b, c), we similarly have We have used Eqs. (A.7) and (A.8). This established the connection between L 011 (a, b, c) and L 000 (a, b, c).
For the completeness, we note that L 000 (a, b, c) have somewhat laconic integral representation In the following, we give an explicit computation of the integral. We first change of variables and introduce z = e iφ ; thus the angular integration over φ from −π to π becomes a closed contour integral of |z| = 1 in counter-clockwise direction.
The integral becomes (for real a, b and with b > 0) To proceed, we have to consider specific relations between a and b, When a > b > 0, the solution z − lies outside of the unit circle |z| = 1. Thus z + is the only pole inside the unit circle. Using the residue theorem, one obtains the result In this case, the situation is more subtle. We can write If = 0, the two roots are located exactly on the unit circle. We are interested in the case of = 0 + . The root with negative imaginary part is slightly shifted inside the unit circle while the solution with positive imaginary part is slightly shifted outside the unit circle. Following the residue theorem, one obtains Amusingly, we can apply this formulae to explicitly derive the identity In this case, z + (z − ) is located outside (inside) of the unit circle. Thus To sum up, we evaluated the auxiliary integral to obtain π −π dφ 2π (B.14) Using Eq. (B.14), one can prove all the formulae involving integrals of three Bessel functions given in Appendix A.
C Analytic structure of the functions I 1,2,3 and G 1,2,3 In this appendix we consider the analytic structure of the functions I 1,2,3 and G 1,2,3 . In particular we concentrate on possible singularities.
• I 1 (p, q, k): First for the readers convenience, we recall the definition The behavior of the function I 1 (p, q, k) depends on the input momenta p, q, k. Let's consider a few cases.
When k × q = 0 or q × p = 0 or simply q = 0, the first term simply vanishes.
When p ⊥ = |q − p|, the angular integration includes the point w ⊥ = 0 1 . Thus it appears that I 1 may have a singularity in this case. We want to show that this is a removable singularity. Consider the limit w ⊥ → 0 of the ratio We thus conclude that this ratio is finite. On the other hand, when p ⊥ = |q − p| and w ⊥ = 0 the numerators in the integrands vanish. This proves that w ⊥ = 0 is a removable singularity of the integrand and the value of the integrand is zero.
When w ⊥ = q ⊥ , one may suspect the first integrand is singular. Consider the piece: We performed Taylor expansion around w ⊥ = q ⊥ for the ratio J 1 (w ⊥ τ )/w ⊥ in the last equality. Thus in the limit w ⊥ → q ⊥ , the integrand becomes We are left with an integral of the form For the case of interest, |b − c| < a < b + c, one can evaluate L 011 (a, b, c): Using this result it is straightforward to compute d db Substituting a = |k−q|, b = q ⊥ , c = k ⊥ into the formula, one can compute Eq. (C.4). to obtain the final expression for the integrand at Recall that from we have implicitly assumed that q × k = 0. Therefore w ⊥ = q ⊥ is a removable pole.
The superficial singularity at w ⊥ = q ⊥ is analysed by The integrand becomes When |a − b| < c < a + b, There are no superficial singularities at w ⊥ = 0 or q ⊥ = 0.

D Comparison with partial results from LCPT
Previously Chirilli, Kovchegov and Wertepny (CKW) derived the single-gluon production amplitude (in general, a gauge-variant object) at leading-order and order-g 3 in Ref. [2]; the full expression for the first saturation correction to the gluon production cross-section (a gauge-invariant object) was not obtained. We thus are forced to compare the leadingorder M 1 and order-g 3 (M 3 ) amplitudes only. However, in contrast to our approach, these calculations were done in the light cone gauge. Naively, this difference in gauges would not allow us to perform one-to-one comparison of the results on the level of the amplitudes as they are gauge-variant objects. Deeper inspection shows that the leading order amplitude is gauge invariant as it is the only contribution to the leading order particle production cross-section (order-by-order gauge invariant-object). At the leading order our calculations are in full agreement with Ref. [2] and previous studies. The first saturation correction contribution to single inclusive gluon production involves leading-order, orderg 3 and order-g 5 amplitudes: M 3 (k)M * 3 (k) + M 1 (k)M * 5 (k) + M * 1 (k)M 5 (k). This is a gauge invariant combination. Since M 5 was not computed in Ref. [2] we cannot perform the comparison with our result. Remarkably we can compare the odd part of M 3 (k): that is M 3 (k) − M 3 (−k). Indeed, the leading odd contribution to the double inclusive gluon production is given by (n (3) (k) − n (3) (−k))(n (3) (p) − n (3) (−p)) , which in the amplitude notation would correspond to M 1 (k)(M 3 (k) − M 3 (−k))M 1 (p)(M 3 (p) − M 3 (−p)) . This is an observable, that is a gauge invariant object. Moreover as we established before, M 1 (k) is gauge invariant, therefore (M 3 (k) − M 3 (−k)) is gauge invariant as well. Now the final subtle point is that the calculations in Ref. [2] were performed not in terms of continuous