aainstitutetext: Department of Physics, North Carolina State University,
Raleigh, NC 27695, USAbbinstitutetext: RIKEN/BNL Research Center, Brookhaven
National Laboratory, Upton, NY 11973

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

Ming Li a,b and Vladimir V. Skokov<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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. Balitsky:2004rr ) and
independently by Chirilli, Kovchegov and Wertepny (see Ref. 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.

##  1 Introduction

In Ref. Li:2021zmf , 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 next-to-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. Balitsky:2004rr ; Chirilli:2015tea  for the single inclusive gluon
production and in Refs. McLerran:2016snu ; Kovchegov:2018jun  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. Chirilli:2015tea , 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 Romatschke:2005pm ; Gelis:2013rba ; Berges:2013eia , 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 Kovchegov:1997pc . The saturation momentum $Q_{s}$
is proportional to $\alpha_{s}\mu$, where $\mu^{2}$ is the hadron color charge
squared per unit transverse area Lappi:2007ku . Therefore, the single
inclusive gluon productions in high energy nuclear collisions can be expressed
as a function of two dimensionless parameters $Q_{s,P}^{2}/k_{\perp}^{2}$ and
$Q_{s,T}^{2}/k_{\perp}^{2}$, involving the saturation scales of the projectile
and the target respectively. In the notations of Ref. Chirilli:2015tea ,

$\frac{dN}{d^{2}\mathbf{k}dy}=\frac{1}{\alpha_{s}}f\left(\frac{Q_{s,P}^{2}}{k_{\perp}^{2}},\frac{Q_{s,T}^{2}}{k_{\perp}^{2}}\right)\,.$
(1)

Here the overall factor of $\alpha_{s}^{-1}$ 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_{s,P}^{2}/k_{\perp}^{2}$:

$f\left(\frac{Q_{s,P}^{2}}{k_{\perp}^{2}},\frac{Q_{s,T}^{2}}{k_{\perp}^{2}}\right)=\sum_{n=1}^{\infty}\left(\frac{Q_{s,P}^{2}}{k_{\perp}^{2}}\right)^{n}f_{n}\left(\frac{Q_{s,T}^{2}}{k_{\perp}^{2}}\right)\,.$
(2)

The coefficients $f_{n}$ are analytically computable and can be non-
perturbatively evaluated to all order in $Q_{s,T}^{2}/k_{\perp}^{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.
Kovchegov:1998bi ; Dumitru:2001ux :

$\frac{dN}{d^{2}\mathbf{k}dy}\approx\frac{1}{\alpha_{s}}\frac{Q_{s,P}^{2}}{k_{\perp}^{2}}f_{1}\left(\frac{Q_{s,T}^{2}}{k_{\perp}^{2}}\right)\,.$
(3)

This truncation is justified as long as one is interested in the momentum
region $k_{\perp}\gg Q_{s,P}$ of the gluon spectrum. As $k_{\perp}$ approaches
$Q_{s,P}$, higher order corrections have to be taken into account. For
$k_{\perp}\leq Q_{s,P}$, an all-order resummation is needed (see e.g. Refs.
Chen:2015wia ; Li:2016eqr , and numerical simulations of Refs. Krasnitz:1999wc
; Krasnitz:2000gz ; Schenke:2015aqa ; Schlichting:2019bvy ). We are mainly
interested in the production of semi-hard gluons with transverse momentum
$k_{\perp}>Q_{s,P}$. Thus our goal is to compute the next-to-leading
correction $f_{2}$ in the expansion (2); the partial result for $f_{2}$ was
presented in Ref. Chirilli:2015tea  where the contributions with amplitudes of
order ${\cal O}(g^{5})$ were not included. In Ref. Chirilli:2015tea , 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. Li:2021zmf . 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 $\rho_{P}$ and the target
$\rho_{T}$. We aim to provide the functional
$\frac{dN}{d^{2}\mathbf{k}dy}\left[\rho_{P},\rho_{T}\right]$.

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 $\frac{dN}{d^{2}\mathbf{k}dy}$; it is often obtained
from $\frac{dN}{d^{2}\mathbf{k}dy}\left[\rho_{P},\rho_{T}\right]$ by using the
Gaussian McLerran-Venugopalan model McLerran:1993ni ; McLerran:1993ka . This
may miss some interesting effects potentially relevant for the observables.
Indeed, we will demonstrate that, in general, the functional
$\frac{dN}{d^{2}\mathbf{k}dy}\left[\rho_{P},\rho_{T}\right]$ is not an even
function of $\rho_{P}$, e.g. it contains a product of three color charge
densities $\rho_{P}(\mathbf{x})\rho_{P}(\mathbf{y})\rho_{P}(\mathbf{z})$. The
terms with odd powers of $\rho_{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 Jeon:2005cf ;
Dumitru:2020fdh ; Dumitru:2021tvw . For instance, the proton model of Ref.
Dumitru:2020fdh  includes a nonzero C-odd part of the correlator of three
color charges. Our resulting functional
$\frac{dN}{d^{2}\mathbf{k}dy}\left[\rho_{P},\rho_{T}\right]$ will explicitly
keep all even and odd powers of $\rho_{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
$\frac{dN}{d^{2}\mathbf{k}dy}\left[\rho_{P},\rho_{T}\right]$ allows us to
evaluate multiparticle production in a very straightforward way Gelis:2008rw ;
Gelis:2008ad ; Gelis:2008sz  – by finding an average of the product of the
functionals, e.g. for two-particle inclusive we simply get

$\frac{d^{2}N}{d^{2}\mathbf{k}_{1}dy_{1}d^{2}\mathbf{k}_{2}dy_{2}}=\left\langle\frac{dN}{d^{2}\mathbf{k}_{1}dy_{1}}\left[\rho_{P},\rho_{T}\right]\frac{dN}{d^{2}\mathbf{k}_{2}dy_{2}}\left[\rho_{P},\rho_{T}\right]\right\rangle_{P,T}$
(4)

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.
Kovner:2010xk ; Kovner:2016jfp  (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\rho^{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.

Figure 1: Schematic diagrams showing the leading order single gluon production
in high energy proton-nucleus collisions. The $\rho^{b}(\mathbf{x})$ is the
valence color charge density of the proton. The colored solid rectangular
represents the highly Lorentz contracted nucleus.

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. Chirilli:2015tea  for the complete 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}\rho^{b_{1}}(\mathbf{x}_{1})\rho^{b_{2}}(\mathbf{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.

Figure 2: Partial diagrams representing the order-$g^{3}$ single gluon
production amplitude. It involves interactions of two color charge densities
$\rho^{b_{1}}(\mathbf{x}_{1})$ and $\rho^{b_{2}}(\mathbf{x}_{2})$.  Figure 3:
Partial diagrams showing the order-$g^{5}$ single gluon production amplitude.
It involves interactions of three color charge densities
$\rho^{b_{1}}(\mathbf{x}_{1})$, $\rho^{b_{2}}(\mathbf{x}_{2})$ and
$\rho^{b_{3}}(\mathbf{x}_{3})$.

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
$g^{5}\rho^{b_{1}}(\mathbf{x}_{1})\rho^{b_{2}}(\mathbf{x}_{2})\rho^{b_{3}}(\mathbf{x}_{3})$.
These diagrams could 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)}+\ldots$. 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)}\sim g\rho_{P}$, $M_{(3)}\sim
g^{3}\rho_{P}^{2}$ and $M_{(5)}\sim g^{5}\rho_{P}^{3}$. The single gluon
production is proportional to the amplitude squared

$|M|^{2}=|M_{(1)}|^{2}+M_{(1)}M_{(3)}^{\ast}+M_{(1)}^{\ast}M_{(3)}+|M_{(3)}|^{2}+M_{(1)}M_{(5)}^{\ast}+M_{(1)}^{\ast}M_{(5)}+\ldots$
(5)

The crossing terms between order-$g$ and order-$g^{3}$ amplitudes are
proportional to odd powers of the color charge density
$M_{(1)}M_{(3)}^{\ast}+M_{(1)}^{\ast}M_{(3)}\sim g^{4}\rho_{P}^{3}$. In
McLerran-Venugopalan model (see also Ref. Kovchegov:1996ty ), $\rho_{P}$
satisfies a Gaussian distribution; therefore $\rho_{P}^{3}$ terms vanish upon
averaging over all possible configurations of color charge density. The first
saturation correction therefore starts from

$\frac{dN}{d^{2}\mathbf{k}dy}\Big{|}_{FSC}\sim|M_{(3)}(\mathbf{k})|^{2}+M_{(1)}(\mathbf{k})M_{(5)}^{\ast}(\mathbf{k})+M_{(1)}^{\ast}(\mathbf{k})M_{(5)}(\mathbf{k}).$
(6)

Both the order-$g^{3}$ amplitude squared and the crossing terms between
order-$g$ and order-$g^{5}$ terms are proportional to $g^{6}\rho_{P}^{4}$. In
a Gaussian model for the color charge density, one has $\langle
g^{4}\rho_{P}^{2}\rangle\sim g^{4}\mu_{P}^{2}\sim Q_{s,P}^{2}$ so that the
first saturation correction is proportional to
$\frac{1}{\alpha_{s}}Q_{s,P}^{4}$. This should be compared with the leading
order term $|M_{(1)}|^{2}\sim\frac{1}{\alpha_{s}}Q_{s,P}^{2}$. This power
counting in terms of the saturation scale has already been formally shown in
Eq. (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 Li:2021zmf . 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
$\frac{dN}{d^{2}\mathbf{k}dy}\left[\rho_{P},\rho_{T}\right]$ 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.

##  2 General Setup and Solutions of Classical Yang-Mills Equations

We briefly recapitulate the general setup of our calculations and the results
obtained in Ref. Li:2021zmf . In the CGC framework, the classical gluon fields
produced in high energy proton-nucleus collisions are determined by the
classical Yang-Mills equations $D_{\mu}F^{\mu\nu}=J^{\nu}$. Here the color
charge current for a right-moving projectile and a left-moving target
$J^{\nu}=\delta^{\nu+}\rho_{P}(x^{-},\mathbf{x})+\delta^{\nu-}\rho_{T}(x^{+},\mathbf{x})$.
In the forward light cone $x^{+}>0$ and $x^{-}>0$, the color charge current
vanishes and the effects of the sources $\rho_{P},\rho_{T}$ are reflected in
the initial conditions. We adopt the procedure described in Refs.
Kovner:1995ja ; Dumitru:2001ux  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,
$A^{+}=A_{-}=x^{+}\alpha(\tau,\mathbf{x})$,
$A^{-}=A_{+}=-x^{-}\alpha(\tau,\mathbf{x})$,
$A^{i}=\alpha^{i}(\tau,\mathbf{x})$ with $\tau=\sqrt{2x^{+}x^{-}}$. In terms
of $\alpha(\tau,\mathbf{x}),\alpha^{i}(\tau,\mathbf{x})$, the classical Yang-
Mills equations read

$\begin{split}&\partial_{\tau}^{2}\alpha+\frac{3}{\tau}\partial_{\tau}\alpha-[D_{i},[D_{i},\alpha]]=0\,,\\\
&-ig[\alpha,\tau\partial_{\tau}\alpha]+[D_{i},\frac{1}{\tau}\partial_{\tau}\alpha_{i}]=0\,,\\\
&\frac{1}{\tau}\partial_{\tau}\alpha_{i}+\partial^{2}_{\tau}\alpha_{i}-ig\tau^{2}[\alpha,[D_{i},\alpha]-D_{j}F_{ji}=0\\\
\end{split}$ (7)

with the initial conditions Kovner:1995ja ; Gyulassy:1997vt :

$\begin{split}&\alpha(\tau=0,\mathbf{x})=\frac{ig}{2}[\alpha_{P}^{i}(\mathbf{x}),\alpha_{T}^{i}(\mathbf{x})],\\\
&\alpha^{i}(\tau=0,\mathbf{x})=\alpha_{P}^{i}(\mathbf{x})+\alpha_{T}^{i}(\mathbf{x}).\\\
\end{split}$ (8)

The $\alpha_{P}^{i}(\mathbf{x})$ and $\alpha_{T}^{i}(\mathbf{x})$ are the
Weizsacker-Williams (WW) gluon fields of the projectile and target,
respectively; they depend on the transverse coordinate and are related to
$\rho_{P(T)}$ through the static Yang-Mills equations
$\partial^{i}\alpha_{P(T)}^{i}(\mathbf{x})=g\rho_{P(T)}$. Since each field
$\alpha_{P(T)}^{i}(\mathbf{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
$\alpha_{P}^{i}(\mathbf{x})=\frac{i}{g}V(\mathbf{x})\partial^{i}V^{\dagger}(\mathbf{x})$,
$\alpha_{T}^{i}(\mathbf{x})=\frac{i}{g}U(\mathbf{x})\partial^{i}U^{\dagger}(\mathbf{x})$.
In general, the WW field of the projectile (target) contains all powers in
$\rho_{P}$ ($\rho_{T}$); for our purpose only a power series expansion in
$\rho_{P}$ is required:

$\alpha_{P}^{i}(\mathbf{x})=\sum_{n=0}^{\infty}g^{2n+1}\alpha_{P}^{{(2n+1)},i}(\tau,\mathbf{x})\,.$
(9)

The first three expansion coefficients are

$\displaystyle\alpha_{P}^{(1),i}$ $\displaystyle=\partial^{i}\phi,$ (10)
$\displaystyle\alpha_{P}^{(3),i}$
$\displaystyle=\frac{1}{2}i\left(\delta^{ij}-\frac{\partial^{i}\partial^{j}}{\partial^{2}}\right)[\phi,\partial^{j}\phi],$
(11) $\displaystyle\alpha_{P}^{(5),i}$
$\displaystyle=\frac{1}{2}\left(\delta^{ij}-\frac{\partial^{i}\partial^{j}}{\partial^{2}}\right)\left([\frac{1}{\partial^{2}}[\phi,\partial^{2}\phi],\partial^{j}\phi]-\frac{1}{3}[\phi,[\phi,\partial^{j}\phi]]\right),$
(12)

where the gluon exchange potential

$\phi(\mathbf{x})=\frac{1}{\partial^{2}}\rho_{P}\equiv\frac{1}{2\pi}\int
d^{2}\mathbf{y}\ln{\left(|\mathbf{x}-\mathbf{y}|m\right)}\rho_{P}(\mathbf{y})\,$
(13)

with $m$ being some appropriate IR cutoff.

The initial value problem can be solved order-by-order in $\rho_{P}$ while
preserving non-perturbative dependence to all orders in $\rho_{T}$. In order
to facilitate this, we consider a subgauge transformation defined by the
condition $\partial^{i}\beta^{i}(\tau=0,\mathbf{x})=0$ to which we will refer
to as the initial time Coulomb gauge, see e.g. Ref. Blaizot:2010kh ;
Berges:2013fga . The fields $\beta$ and $\beta^{i}$ are non-linearly related
to $\alpha$ and $\alpha_{i}$ through a gauge transformation

$\begin{split}&\beta^{i}=\Omega^{+}\alpha^{i}\Omega+\frac{i}{g}\Omega^{+}\partial^{i}\Omega,\\\
&\beta=\Omega^{+}\alpha\Omega.\\\ \end{split}$ (14)

The $\Omega=\Omega(\mathbf{x})$ is fixed to satisfy
$\partial^{i}\beta^{i}(\tau=0,\mathbf{x})=0$. A closed form of $\Omega$ is not
known but a perturbative expansion in coupling constant is obtained in Ref.
Li:2021zmf . 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\rho_{T}$). Due to the structure of the
classical Yang-Mills equations, only odd powers of $g$ are non-vanishing, that
is

$\begin{split}&\beta(\tau,\mathbf{x})=\sum_{n=0}^{\infty}g^{2n+1}\beta^{{(2n+1)}}(\tau,\mathbf{x}),\\\
&\beta_{i}(\tau,\mathbf{x})=\sum_{n=0}^{\infty}g^{2n+1}\beta_{i}^{{(2n+1)}}(\tau,\mathbf{x}).\\\
\end{split}$ (15)

The order-$g$ and order-$g^{3}$ solutions in momentum space are

$\begin{split}&\beta^{(1)}(\tau,\mathbf{k})=b_{\eta}(\mathbf{k})\frac{J_{1}(k_{\perp}\tau)}{k_{\perp}\tau},\\\
&\beta^{(1)}_{i}(\tau,\mathbf{k})=\frac{-i\epsilon^{il}\mathbf{k}_{l}}{k_{\perp}^{2}}b_{\perp}(\mathbf{k})J_{0}(k_{\perp}\tau)\\\
\end{split}$ (16)

and

$\begin{split}\beta^{(3)}(\tau,\mathbf{q})=&2\beta^{(3)}(\tau=0,\mathbf{q})\frac{J_{1}(q_{\perp}\tau)}{q_{\perp}\tau}+i\int\frac{d^{2}\mathbf{p}}{(2\pi)}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\Big{[}b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\Big{]}\\\
&\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})+w_{\perp}^{2}}{q_{\perp}^{2}-w_{\perp}^{2}}\left(\frac{J_{1}(w_{\perp}\tau)}{w_{\perp}\tau}-\frac{J_{1}(q_{\perp}\tau)}{q_{\perp}\tau}\right),\\\
\end{split}$ (17)
$\begin{split}\beta^{(3)}_{\perp}(\tau,\mathbf{q})=&\beta_{\perp}^{(3)}(\tau=0,\mathbf{q})J_{0}(q_{\perp}\tau)+\frac{i}{q_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p^{2}_{\perp}|\mathbf{q}-\mathbf{p}|^{2}}\Big{(}\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\Big{[}b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\Big{]}\\\
&+(\mathbf{p}\cdot\mathbf{q}-p_{\perp}^{2}-q_{\perp}^{2})\Big{[}b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})\Big{]}\Big{)}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}-q_{\perp}^{2}}(J_{0}(w_{\perp}\tau)-J_{0}(q_{\perp}\tau))\\\
&+\frac{i}{q_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{q}\times\mathbf{p}}{2p^{2}_{\perp}|\mathbf{q}-\mathbf{p}|^{2}}\Big{[}b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\Big{]}\Big{(}J_{0}(p_{\perp}\tau)J_{0}(|\mathbf{q}-\mathbf{p}|\tau)-J_{0}(q_{\perp}\tau)\Big{)},\end{split}$
(18)
$\begin{split}\Lambda^{(3)}(\tau,\mathbf{q})=&-\frac{i}{q_{\perp}^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})}{4p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\Big{(}2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\left[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})\right]\\\
\
&-(p_{\perp}^{2}+|\mathbf{q}-\mathbf{p}|^{2})\left[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\right]\Big{)}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w^{2}_{\perp}}(1-J_{0}(w_{\perp}\tau))\\\
&-\frac{i}{q_{\perp}^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})}{4p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\left[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\right]\Big{(}1-J_{0}(p_{\perp}\tau)J_{0}(|\mathbf{q}-\mathbf{p}|\tau)\Big{)}.\end{split}$
(19)

Here
$w^{2}_{\perp}=p^{2}_{\perp}+|\mathbf{q}-\mathbf{p}|^{2}-2p_{\perp}|\mathbf{q}-\mathbf{p}|\cos{\phi}$.
The decompositions of transverse field in momentum space is

$\beta_{i}^{(3)}(\tau,\mathbf{q})=\frac{-i\epsilon^{ij}\mathbf{q}_{j}}{q_{\perp}}\beta_{\perp}^{(3)}(\tau,\mathbf{q})-i\mathbf{q}_{i}\Lambda^{(3)}(\tau,\mathbf{q}).$
(20)

The above solutions are expressed in terms of order-$g$ initial conditions
$b_{\eta}(\mathbf{k})$, $b_{\perp}(\mathbf{k})$ and order-$g^{3}$ initial
conditions $\beta^{(3)}(\tau=0,\mathbf{q})$,
$\beta_{\perp}^{(3)}(\tau=0,\mathbf{q})$. To present expressions for the
initial conditions, we introduce an auxiliary quantity which will repeatedly
appear in our calculations:

$\zeta^{i}(\mathbf{x})=U^{\dagger}(\mathbf{x})\alpha_{P}^{i}(\mathbf{x})U(\mathbf{x})=\alpha_{P}^{a,i}(\mathbf{x})U^{ac}(\mathbf{x})T^{c}\,.$
(21)

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 $\zeta^{i}(\mathbf{x})$ in terms of coupling constant $g$
with the coefficients

$\zeta^{i}_{(n)}(\mathbf{x})=U^{\dagger}(\mathbf{x})\alpha_{P}^{(n),i}(\mathbf{x})U(\mathbf{x})\,.$
(22)

For example, the explicit expression for $\zeta^{i}(\mathbf{x})$ at order-$g$
is

$\zeta^{i}_{(1)}(\mathbf{k})=\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{i\mathbf{p}_{i}}{p_{\perp}^{2}}\rho_{P}^{a}(\mathbf{p})U^{ac}(\mathbf{k}-\mathbf{p})T^{c}\,.$
(23)

By analogy with $\zeta^{i}(\mathbf{x})$, we also introduce the eikonally
rotated projectile color charge density

$\rho_{R}(\mathbf{x})=U^{\dagger}(\mathbf{x})\partial^{i}\alpha_{P}^{(1),i}U(\mathbf{x})=g\,U^{\dagger}(\mathbf{x})\rho_{P}(\mathbf{x})U(\mathbf{x})\,.$
(24)

In terms of these auxiliary quantities, we have order-$g$ initial conditions

$\begin{split}&b_{\eta}(\mathbf{k})=2\beta^{(1)}(\tau=0,\mathbf{k})=-i\mathbf{k}^{i}\zeta^{(1),i}(\mathbf{k})-g\rho_{R}(\mathbf{k}),\\\
&b_{\perp}(\mathbf{k})=-i\epsilon_{ij}\mathbf{k}_{i}\beta_{j}^{(1)}(\tau=0,\mathbf{k})=-i\epsilon^{ij}\mathbf{k}^{i}\zeta^{(1),j}(\mathbf{k}).\\\
\end{split}$ (25)

In order to make the further notation uniform we also introduce

$b_{\|}(\mathbf{k})=\frac{i\mathbf{k}^{i}}{k_{\perp}^{2}}\zeta^{(1),i}(\mathbf{k}).$
(26)

In other words we can expand the two dimensional vector $\zeta^{(1),i}$ into
the longitudinal and transverse parts

$\zeta^{i}_{(1)}(\mathbf{k})=\frac{-i\epsilon^{ij}\mathbf{k}^{j}}{k_{\perp}^{2}}b_{\perp}(\mathbf{k})-i\mathbf{k}^{i}b_{\|}(\mathbf{k}).$
(27)

The order-$g^{3}$ initial conditions are

$\begin{split}\beta_{\perp}^{(3)}(\tau=0,\mathbf{k})=&\frac{i\epsilon_{ih}\mathbf{k}_{h}}{k_{\perp}}\zeta_{(3)}^{i}(\mathbf{k})-i\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot\mathbf{p}}{k_{\perp}p_{\perp}^{2}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),b_{\perp}(\mathbf{p})\right]\\\
&-i\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{p}\times\mathbf{k}}{2k_{\perp}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),b_{\|}(\mathbf{p})\right],\end{split}$
(28)
$2\beta_{(3)}(\tau=0,\mathbf{k})=-i\mathbf{k}^{i}\zeta_{(3)}^{i}(\mathbf{k})-i\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),b_{\eta}(\mathbf{p})\right].$
(29)

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

$\begin{split}\beta_{\perp}^{(5)}(\tau=0,\mathbf{k})=&\frac{i\epsilon_{lh}\mathbf{k}_{h}}{k_{\perp}}\zeta^{l}_{(5)}(\mathbf{k})+\frac{\epsilon^{lh}\mathbf{k}^{h}}{k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),\zeta^{l}_{(3)}(\mathbf{p})\right]\\\
&+\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{k_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}^{l}}{p_{\perp}^{2}}\left[\zeta^{l}_{(3)}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\right]\\\
&-\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{k_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}\times\mathbf{q}}{p_{\perp}^{2}q_{\perp}^{2}}\left[\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\perp}(\mathbf{q})\right],b_{\perp}(\mathbf{k}-\mathbf{p})\right]\\\
&-\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{k_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}\cdot\mathbf{q}}{2p_{\perp}^{2}}\left[\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\|}(\mathbf{q})\right],b_{\perp}(\mathbf{k}-\mathbf{p})\right]\\\
&-\frac{1}{2}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot\mathbf{q}}{k_{\perp}q_{\perp}^{2}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\perp}(\mathbf{q})\right]\right]\\\
&-\frac{1}{3}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{q}\times\mathbf{k}}{k_{\perp}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\|}(\mathbf{q})\right]\right].\\\
\end{split}$ (30) $\begin{split}&2\beta_{(5)}(\tau=0,\mathbf{k})\\\
=&-i\mathbf{k}^{i}\zeta_{(5)}^{i}(\mathbf{k})-\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),\mathbf{p}^{i}\zeta_{(3)}^{i}(\mathbf{p})\right]+\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{p}^{l}}{p_{\perp}^{2}}\left[\zeta_{(3)}^{l}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\right]\\\
&-\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{p}\times\mathbf{q}}{p_{\perp}^{2}q_{\perp}^{2}}\left[\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\perp}(\mathbf{q})\right],b_{\eta}(\mathbf{k}-\mathbf{p})\right]\\\
&-\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{p}\cdot\mathbf{q}}{2p_{\perp}^{2}}\left[\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\|}(\mathbf{q})\right],b_{\eta}(\mathbf{k}-\mathbf{p})\right]\\\
&-\frac{1}{2}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\left[b_{\|}(\mathbf{k}-\mathbf{p}),\left[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\eta}(\mathbf{q})\right]\right].\\\
\end{split}$ (31)

##  3 Defining the Semi-Hard Gluon Productions

In this section we conform the conventional four-dimensional Minkowski LSZ
procedure (see the textbook Itzykson:1980rh ) to the problem at hand, that is
2+1 Milne coordinate system with the initial conditions defined at $\tau=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
$\tau=0$.

###  3.1 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 $\beta(\tau,\mathbf{x})$ and
$\beta_{\perp}(\tau,\mathbf{x})$. The remaining field
$\Lambda(\tau,\mathbf{x})$ is non-dynamical and does not contribute to the
particle spectrum. The free-field equations for
$\tilde{\beta}(\tau,\mathbf{x})\equiv\tau\beta(\tau,\mathbf{x})$ and
$\beta_{\perp}(\tau,\mathbf{x})$ are

$\begin{split}&\tau^{2}\partial_{\tau}^{2}\tilde{\beta}+\tau\partial_{\tau}\tilde{\beta}-\tilde{\beta}-\tau^{2}\partial^{2}_{i}\tilde{\beta}=0\,,\\\
&\tau^{2}\partial^{2}_{\tau}\beta_{\perp}+\tau\partial_{\tau}\beta_{\perp}-\tau^{2}\partial^{2}\beta_{\perp}=0.\\\
\end{split}$ (32)

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 $\tau$. This fixes the form of the functions – they are
Hankel functions of first and second kind (see Ref. abramowitz+stegun  for
mathematical definition); the first (zeroes) order Hankel function solves the
first (second) equation in Eqs. (32). 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

$f(k_{\perp}\tau)=\sqrt{\frac{\pi}{4}}H_{1}^{(1)}(k_{\perp}\tau),\qquad
f^{\ast}(k_{\perp}\tau)=\sqrt{\frac{\pi}{4}}H_{1}^{(2)}(k_{\perp}\tau)$ (33)

and

$g(k_{\perp}\tau)=\sqrt{\frac{\pi}{4}}H_{0}^{(1)}(k_{\perp}\tau),\qquad
g^{\ast}(k_{\perp}\tau)=\sqrt{\frac{\pi}{4}}H_{0}^{(2)}(k_{\perp}\tau).$ (34)

They satisfy the Wronskian identities

$\begin{split}&f(k_{\perp}\tau)\partial_{\tau}f^{\ast}(k_{\perp}\tau)-\partial_{\tau}f(k_{\perp}\tau)f^{\ast}(k_{\perp}\tau)=\frac{-i}{\tau},\\\
&g(k_{\perp}\tau)\partial_{\tau}g^{\ast}(k_{\perp}\tau)-\partial_{\tau}g(k_{\perp}\tau)g^{\ast}(k_{\perp}\tau)=\frac{-i}{\tau}.\\\
\end{split}$ (35)

Next we construct the mode expansions for both dynamical fields in terms of
creation and annihilation operators

$\begin{split}\tilde{\beta}^{a}(\tau,\mathbf{x})&=\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\left(\hat{a}^{a}_{\mathbf{k}}f^{\ast}(k_{\perp}\tau)e^{i\mathbf{k}\cdot\mathbf{x}}+\hat{a}_{\mathbf{k}}^{a\dagger}f(k_{\perp}\tau)e^{-i\mathbf{k}\cdot\mathbf{x}}\right),\\\
\beta_{\perp}^{a}(\tau,\mathbf{x})&=\int\frac{d^{2}\mathbf{k}}{(2\pi)^{2}}\left(\hat{c}_{\mathbf{k}}^{a}g^{\ast}(k_{\perp}\tau)e^{i\mathbf{k}\cdot\mathbf{x}}+\hat{c}_{\mathbf{k}}^{a\dagger}g(k_{\perp}\tau)e^{-i\mathbf{k}\cdot\mathbf{x}}\right)\,.\end{split}$
(36)

One can check that the normalization of the mode functions leads to the
consistent commutation relations, that is

$\begin{split}&[\tilde{\beta}^{a}(\tau,\mathbf{x}),\tau\partial_{\tau}\tilde{\beta}^{b}(\tau,\mathbf{y})]=i\delta^{ab}\delta^{(2)}(\mathbf{x}-\mathbf{y}),\\\
&[\beta_{\perp}^{a}(\tau,\mathbf{x}),\tau\partial_{\tau}\beta_{\perp}^{b}(\tau,\mathbf{y})]=i\delta^{ab}\delta^{(2)}(\mathbf{x}-\mathbf{y})\\\
\end{split}$ (37)

for the fields and

$\begin{split}&\left[\hat{a}^{a}_{\mathbf{k}},\hat{a}^{b\dagger}_{\mathbf{p}}\right]=(2\pi)^{2}\delta^{(2)}(\mathbf{k}-\mathbf{p})\delta^{ab},\\\
&\left[\hat{c}^{a}_{\mathbf{k}},\hat{c}^{b\dagger}_{\mathbf{p}}\right]=(2\pi)^{2}\delta^{(2)}(\mathbf{k}-\mathbf{p})\delta^{ab}\\\
\end{split}$ (38)

for the creation and annihilation operators. The mode expansion (36) can be
inverted in order to express the creation operators in terms of the fields

$\begin{split}\hat{a}_{\mathbf{p}}^{a\dagger}=&-i\tau\sqrt{\frac{\pi}{4}}\left(H_{1}^{(2)}(p_{\perp}\tau)\overleftrightarrow{\partial_{\tau}}\tilde{\beta}^{a}(\tau,\mathbf{p})\right),\\\
\hat{c}_{\mathbf{p}}^{a\dagger}=&-i\tau\sqrt{\frac{\pi}{4}}\left(H_{0}^{(2)}(p_{\perp}\tau)\overleftrightarrow{\partial_{\tau}}\beta^{a}_{\perp}(\tau,\mathbf{p})\right).\\\
\end{split}$ (39)

Here
$f_{1}(x)\overleftrightarrow{\partial_{x}}f_{2}(x)=f_{1}(x)\partial_{x}f_{2}(x)-\partial_{x}f_{1}(x)f_{2}(x)$.
Similar to Minkovski spacetime, we can start from the trivial identity for the
creation operator for the asymptotic out state $\tau\rightarrow\infty$:

$\begin{split}\hat{a}_{\mathbf{p}}^{a\dagger}(+\infty)=&\hat{a}_{\mathbf{p}}^{a\dagger}(\tau=0)+\int_{0}^{\infty}d\tau\partial_{\tau}\hat{a}_{\mathbf{p}}^{a\dagger}(\tau)\\\
=&\hat{a}_{\mathbf{p}}^{a\dagger}(\tau=0)-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{1}^{(2)}(p_{\perp}\tau)S_{\eta}^{a}(\tau,\mathbf{p}),\end{split}$
(40)
$\begin{split}\hat{c}_{\mathbf{p}}^{a\dagger}(+\infty)=&\hat{c}_{\mathbf{p}}^{a\dagger}(\tau=0)+\int_{0}^{\infty}d\tau\,\partial_{\tau}\hat{c}_{\mathbf{p}}^{a\dagger}(\tau)\\\
=&\hat{c}_{\mathbf{p}}^{a\dagger}(\tau=0)-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{0}^{(2)}(p_{\perp}\tau)S_{\perp}^{a}(\tau,\mathbf{p}).\\\
\end{split}$ (41)

Here we used the fact that $H_{1}^{(2)}$ ($H_{0}^{(2)}$) satisfies the first
(second) free field equation in Eqs. (32). Additionally we took into account
that the fields $\tilde{\beta}(\tau,\mathbf{p})$ and
$\beta_{\perp}(\tau,\mathbf{p})$ are the (formal) solutions of the non-linear
Yang-Mills equations with the source terms $S_{\eta}(\tau,\mathbf{p})$ and
$S_{\perp}(\tau,\mathbf{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 $\tau=0$, while the second due to the “bulk”
evolution in the forward light cone.

We can also compute $\hat{a}_{p}^{\dagger}(\tau=0)$ and
$\hat{c}_{p}^{\dagger}(\tau=0)$:

$\begin{split}\hat{a}_{\mathbf{p}}^{\dagger}(\tau=0)=&-i\sqrt{\frac{\pi}{4}}\tau\left(H_{1}^{(2)}(p_{\perp}\tau)\partial_{\tau}\tilde{\beta}(\tau,\mathbf{p})-\partial_{\tau}H_{1}^{(2)}(p_{\perp}\tau)\tilde{\beta}(\tau,\mathbf{p})\right)\Big{|}_{\tau\rightarrow
0}\\\ =&C_{1}\sqrt{\frac{1}{\pi}},\end{split}$ (42)
$\begin{split}\hat{c}_{\mathbf{p}}^{\dagger}(\tau=0)=&-i\sqrt{\frac{\pi}{4}}\tau\left(H_{0}^{(2)}(p_{\perp}\tau)\partial_{\tau}\beta_{\perp}(\tau,\mathbf{p})-\partial_{\tau}H_{0}^{(2)}(p_{\perp}\tau)\beta_{\perp}(\tau,\mathbf{p})\right)\Big{|}_{\tau\rightarrow
0}\\\ =&D_{1}\sqrt{\frac{1}{\pi}}.\\\ \end{split}$ (43)

Where the constants $C_{1}$ and $D_{1}$ can be found using the initial
conditions $\beta(\tau=0,\mathbf{p})=\frac{p_{\perp}}{2}C_{1}$ and
$\beta_{\perp}(\tau=0,\mathbf{p})=D_{1}$.

Following the notations of Ref. McLerran:2016snu  we define the surface
contributions $\mathfrak{S}_{\eta}$ and $\mathfrak{S}_{\perp}$ by

$\begin{split}&\mathfrak{S}_{\eta}(\mathbf{p})=\hat{a}_{\mathbf{p}}^{\dagger}(\tau=0)=\frac{1}{\sqrt{\pi}p_{\perp}}2\beta(\tau=0,\mathbf{p}),\\\
&\mathfrak{S}_{\perp}(\mathbf{p})=\hat{c}_{\mathbf{p}}^{\dagger}(\tau=0)=\frac{1}{\sqrt{\pi}}\beta_{\perp}(\tau=0,\mathbf{p}).\\\
\end{split}$ (44)

and the bulk terms $\mathfrak{B}_{\eta}$ and $\mathfrak{B}_{\perp}$ by

$\begin{split}&\mathfrak{B}_{\eta}(\mathbf{p})=-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{1}^{(2)}(p_{\perp}\tau)S_{\eta}(\tau,\mathbf{p}),\\\
&\mathfrak{B}_{\perp}(\mathbf{p})=-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{0}^{(2)}(p_{\perp}\tau)S_{\perp}(\tau,\mathbf{p}).\\\
\end{split}$ (45)

Thus the LSZ reduction formula can be written in the compact form

$\begin{split}&\hat{a}_{\mathbf{p}}^{\dagger}=\mathfrak{S}_{\eta}(\mathbf{p})+\mathfrak{B}_{\eta}(\mathbf{p}),\\\
&\hat{c}_{\mathbf{p}}^{\dagger}=\mathfrak{S}_{\perp}(\mathbf{p})+\mathfrak{B}_{\perp}(\mathbf{p}).\\\
\end{split}$ (46)

As the reader will see below, differentiating between the surface (fully
defined by the initial conditions at $\tau=0$) and the bulk terms (determined
by the final state interactions) helps to pin-point the importance of the
final state interactions.

###  3.2 Single inclusive gluon spectrum

Using the creation and annihilation we define the single inclusive gluon
spectrum (for the fixed valence charge densities) as

$\begin{split}\frac{dN}{d^{2}\mathbf{k}dy}=&\frac{1}{(2\pi)^{2}}\Big{(}\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}+\hat{c}_{\mathbf{k}}^{\dagger}\hat{c}_{\mathbf{k}}\Big{)}\\\
=&\frac{1}{(2\pi)^{2}}\sum_{\gamma=\eta,\perp}\Big{(}|\mathfrak{S}_{\gamma}(\mathbf{k})|^{2}+\mathfrak{S}_{\gamma}(\mathbf{k})\mathfrak{B}_{\gamma}^{\ast}(\mathbf{k})+\mathfrak{S}^{\ast}_{\gamma}(\mathbf{k})\mathfrak{B}_{\gamma}(\mathbf{k})+|\mathfrak{B}_{\gamma}(\mathbf{k})|^{2}\Big{)}\,.\\\
\end{split}$ (47)

Here $\gamma$ denotes two gluon polarizations $\gamma=\eta,\perp$; we sum over
the polarizations as we are only interested in the polarization-blind single
inclusive gluon production.

It is instructive to analyze the $\mathbf{k}\leftrightarrow-\mathbf{k}$
symmetry property of each term in $dN/d^{2}\mathbf{k}dy$ 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 $\mathbf{x}$; consequently the
surface terms have the following symmetry in the momentum space

$\begin{split}&\mathfrak{S}_{\eta}^{a,\ast}(\mathbf{k})=\mathfrak{S}_{\eta}^{a}(-\mathbf{k}),\\\
&\mathfrak{S}_{\perp}^{a,\ast}(\mathbf{k})=\mathfrak{S}_{\perp}^{a}(-\mathbf{k}).\\\
\end{split}$ (48)

The bulk terms require a more detailed analysis. Superficially, one can argue
that since the source terms $S^{a}_{\eta}(\tau,\mathbf{x})$,
$S^{a}_{\perp}(\tau,\mathbf{x})$ are real functions of $\mathbf{x}$ and
consequently
$S^{a,\ast}_{\eta(\perp)}(\tau,\mathbf{k})=S^{a}_{\eta(\perp)}(\tau,-\mathbf{k})$,
we should expect that the entire bulk contribution has to be symmetric under
$\mathbf{k}\leftrightarrow-\mathbf{k}$. This argument however misses the fact
that the Hankel functions in the definition (45) are complex valued. Taking
this into account and explicitly separating the Hankel functions into the real
and imaginary parts we can write

$\begin{split}\mathfrak{B}_{\eta}^{a}(\mathbf{k})=&-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}\frac{d\tau}{\tau}J_{1}(k_{\perp}\tau)S_{\eta}^{a}(\tau,\mathbf{k})-\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}\frac{d\tau}{\tau}Y_{1}(k_{\perp}\tau)S_{\eta}^{a}(\tau,\mathbf{k})\\\
=&\widetilde{\mathfrak{B}}_{\eta}^{a}(\mathbf{k})+\overline{\mathfrak{B}}_{\eta}^{a}(\mathbf{k}).\\\
\end{split}$ (49)

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 $\pi/2$, see a related discussion in Ref. Kovchegov:2018jun . The
newly introduced functions satisfy

$\begin{split}&\widetilde{\mathfrak{B}}_{\eta}^{a,\ast}(\mathbf{k})=-\widetilde{\mathfrak{B}}_{\eta}^{a}(-\mathbf{k}),\\\
&\overline{\mathfrak{B}}_{\eta}^{a,\ast}(\mathbf{k})=\overline{\mathfrak{B}}_{\eta}^{a}(-\mathbf{k}).\end{split}$
(50)

Similar separation can be defined for the transverse component of the bulk
term

$\begin{split}\mathfrak{B}_{\perp}^{a}(\mathbf{k})=&-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}\frac{d\tau}{\tau}J_{0}(k_{\perp}\tau)S_{\perp}^{a}(\tau,\mathbf{k})-\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}\frac{d\tau}{\tau}Y_{0}(k_{\perp}\tau)S_{\perp}^{a}(\tau,\mathbf{k})\\\
=&\widetilde{\mathfrak{B}}_{\perp}^{a}(\mathbf{k})+\overline{\mathfrak{B}}_{\perp}^{a}(\mathbf{k})\\\
\end{split}$ (51)

with the properties

$\begin{split}&\widetilde{\mathfrak{B}}_{\perp}^{a,\ast}(\mathbf{k})=-\widetilde{\mathfrak{B}}_{\perp}^{a}(-\mathbf{k}),\\\
&\overline{\mathfrak{B}}_{\perp}^{a,\ast}(\mathbf{k})=\overline{\mathfrak{B}}_{\perp}^{a}(-\mathbf{k}).\end{split}$
(52)

We are now in position to discuss the symmetry properties of each term in Eq.
(47). The pure surface term is symmetric with respect to
$\mathbf{k}\leftrightarrow-\mathbf{k}$.

$|\mathfrak{S}_{\gamma}(\mathbf{k})|^{2}=\mathfrak{S}^{a,\ast}_{\gamma}(\mathbf{k})\mathfrak{S}^{a}_{\gamma}(\mathbf{k})=\mathfrak{S}^{a}_{\gamma}(-\mathbf{k})\mathfrak{S}^{a}_{\gamma}(\mathbf{k}).$
(53)

The interference terms between the surface and bulk contributions have two
contributions with different properties under
$\mathbf{k}\leftrightarrow-\mathbf{k}$

$\begin{split}&\mathfrak{S}_{\gamma}(\mathbf{k})\mathfrak{B}_{\gamma}^{\ast}(\mathbf{k})+\mathfrak{S}^{\ast}_{\gamma}(\mathbf{k})\mathfrak{B}_{\gamma}(\mathbf{k})\\\
=&\left(-\mathfrak{S}^{a}_{\gamma}(\mathbf{k})\widetilde{\mathfrak{B}}_{\gamma}^{a}(-\mathbf{k})+\mathfrak{S}^{a}_{\gamma}(-\mathbf{k})\widetilde{\mathfrak{B}}^{a}_{\gamma}(\mathbf{k})\right)+\left(\mathfrak{S}^{a}_{\gamma}(\mathbf{k})\overline{\mathfrak{B}}^{a}_{\gamma}(-\mathbf{k}))+\mathfrak{S}^{a}_{\gamma}(-\mathbf{k})\overline{\mathfrak{B}}^{a}_{\gamma}(\mathbf{k}))\right).\\\
\end{split}$ (54)

We grouped the terms according to the symmetry: the first (second) piece
involving $\widetilde{\mathfrak{B}}_{\gamma}^{a}$
($\overline{\mathfrak{B}}_{\gamma}^{a}$) is antisymmetric (symmetric) with
respect to $\mathbf{k}\leftrightarrow-\mathbf{k}$. As was identified in Ref.
McLerran:2016snu , 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

$\begin{split}&|\mathfrak{B}_{\gamma}(\mathbf{k})|^{2}=\mathfrak{B}^{a,\ast}_{\gamma}(\mathbf{p})\mathfrak{B}^{a}_{\gamma}(\mathbf{k})\\\
=&\left(-\widetilde{\mathfrak{B}}^{a}_{\gamma}(-\mathbf{k})\widetilde{\mathfrak{B}}^{a}_{\gamma}(\mathbf{k})+\overline{\mathfrak{B}}^{a}_{\gamma}(-\mathbf{k})\overline{\mathfrak{B}}^{a}_{\gamma}(\mathbf{k})\right)+\left(-\widetilde{\mathfrak{B}}^{a}_{\gamma}(-\mathbf{k})\overline{\mathfrak{B}}^{a}_{\gamma}(\mathbf{k})+\widetilde{\mathfrak{B}}^{a}_{\gamma}(\mathbf{k})\overline{\mathfrak{B}}^{a}_{\gamma}(-\mathbf{k})\right).\end{split}$
(55)

The contribution in the first bracket is symmetric while in the second bracket
is antisymmetric under $\mathbf{k}\leftrightarrow-\mathbf{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 $\frac{dN}{d^{2}\mathbf{k}dy}$ is real.

In practice, our calculations will be done for power series expansion for the
single gluon distribution. We have

$\begin{split}&\mathfrak{S}_{\gamma}=g\mathfrak{S}_{\gamma}^{(1)}+g^{3}\mathfrak{S}_{\gamma}^{(3)}+g^{5}\mathfrak{S}_{\gamma}^{(5)}+\ldots\\\
&\mathfrak{B}_{\gamma}=g^{3}\mathfrak{B}_{\gamma}^{(3)}+g^{5}\mathfrak{B}_{\gamma}^{(5)}+\ldots\\\
\end{split}$ (56)

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}$:

$\begin{split}\frac{dN}{d^{2}\mathbf{k}dy}=&\frac{1}{(2\pi)^{2}}\sum_{\gamma=\eta,\perp}\Bigg{(}g^{2}\left[|\mathfrak{S}_{\gamma}^{(1)}(\mathbf{k})|^{2}\right]\\\
&+g^{4}\Big{[}\mathfrak{S}_{\gamma}^{(1)\ast}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(3)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)}(\mathbf{k})\right)+\mathfrak{S}_{\gamma}^{(1)}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(3)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)\ast}(\mathbf{k})\right)\Big{]}\\\
&+g^{6}\Big{[}|\mathfrak{S}_{\gamma}^{(3)}(\mathbf{k})|^{2}+|\mathfrak{B}_{\gamma}^{(3)}(\mathbf{k})|^{2}+\mathfrak{S}_{\gamma}^{(1)\ast}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(5)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(5)}(\mathbf{k})\right)\\\
&+\mathfrak{S}_{\gamma}^{(1)}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(5)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(5)\ast}(\mathbf{k})\right)+\mathfrak{S}_{\gamma}^{(3)}(\mathbf{k})\mathfrak{B}_{\gamma}^{(3)\ast}(\mathbf{k})+\mathfrak{S}_{\gamma}^{(3)\ast}(\mathbf{k})\mathfrak{B}_{\gamma}^{(3)}(\mathbf{k})\Big{]}\Bigg{)}\\\
\equiv&\frac{1}{(2\pi)^{2}}\Big{(}g^{2}n^{(2)}(\mathbf{k})+g^{4}n^{(4)}(\mathbf{k})+g^{6}n^{(6)}(\mathbf{k})\Big{)}.\end{split}$
(57)

In the following we calculate each order $n^{(i)}(\mathbf{k})$ separately.
From our previous analysis, it is obvious that, for a fixed configuration of
the valence charges, the leading order $n^{(2)}(\mathbf{k})$ is even under
$\mathbf{k}\leftrightarrow-\mathbf{k}$ while all higher order terms have both
even and odd contributions. Note that $g^{4}$ terms is cubic in $\rho_{P}$ and
thus, at least in the Gaussian MV model, does not contribute to the
configuration/event averaged single inclusive gluon production.

##  4 The Gluon Production Amplitude

We now proceed with assembling all the pieces together in order to find the
functional $\frac{dN}{d^{2}\mathbf{k}dy}[\rho_{P},\rho_{T}]$. First we
explicitly write down the surface and the bulk terms.

###  4.1 The surface terms

The surface terms are defined by the initial condition on the light cone
$\tau=0$; thus to find them one does not require solving the classical Yang-
Mills equations in the forward light cone $\tau>0$. Starting from the
definitions

$\begin{split}&\mathfrak{S}_{\eta}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}2\beta(\tau=0,\mathbf{k}),\\\
&\mathfrak{S}_{\perp}(\mathbf{k})=\frac{1}{\sqrt{\pi}}\beta_{\perp}(\tau=0,\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}i\epsilon_{ih}\mathbf{k}_{h}\beta_{i}(\tau=0,\mathbf{k})\\\
\end{split}$ (58)

we obtain that the order-$g$ surface terms are simply

$\begin{split}&\mathfrak{S}_{\eta}^{(1)}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}b_{\eta}(\mathbf{k}),\\\
&\mathfrak{S}_{\perp}^{(1)}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}b_{\perp}(\mathbf{k}).\\\
\end{split}$ (59)

The order-$g^{3}$ and order-$g^{5}$ contribution to the surface terms can be
also extracted straightforwardly

$\mathfrak{S}_{\eta}^{(3)}=\frac{-i\mathbf{k}^{i}}{\sqrt{\pi}k_{\perp}}\zeta_{(3)}^{i}(\mathbf{k})+\frac{ig}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)}\left[b_{\eta}(\mathbf{p}),b_{\|}(\mathbf{k}-\mathbf{p})\right],$
(60)
$\begin{split}\mathfrak{S}_{\perp}^{(3)}=&\frac{i\epsilon^{ih}\mathbf{k}^{h}}{\sqrt{\pi}k_{\perp}}\zeta_{(3)}^{i}(\mathbf{k})+\frac{ig}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot\mathbf{p}}{p_{\perp}^{2}}[b_{\perp}(\mathbf{p}),b_{\|}(\mathbf{k}-\mathbf{p})]\\\
&\qquad-\frac{1}{2}\frac{ig}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}(\mathbf{k}\times\mathbf{p})[b_{\|}(\mathbf{p}),b_{\|}(\mathbf{k}-\mathbf{p})]\\\
\end{split}$ (61)

and

$\begin{split}\mathfrak{S}^{(5)}_{\perp}(\mathbf{k})=&\frac{i\epsilon_{lh}\mathbf{k}_{h}}{\sqrt{\pi}k_{\perp}}\zeta^{l}_{(5)}(\mathbf{k})+\frac{\epsilon^{lh}\mathbf{k}^{h}}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}[b_{\|}(\mathbf{k}-\mathbf{p}),\zeta^{l}_{(3)}(\mathbf{p})]\\\
&+\frac{1}{\sqrt{\pi}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{k_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}^{l}}{p_{\perp}^{2}}\left[\zeta^{l}_{(3)}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\right]\\\
&-\frac{1}{\sqrt{\pi}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{k_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}\times\mathbf{q}}{p_{\perp}^{2}q_{\perp}^{2}}[[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\perp}(\mathbf{q})],b_{\perp}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{\sqrt{\pi}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{k_{\perp}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}\cdot\mathbf{q}}{2p_{\perp}^{2}}[[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\|}(\mathbf{q})],b_{\perp}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{2\sqrt{\pi}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot\mathbf{q}}{k_{\perp}q_{\perp}^{2}}[b_{\|}(\mathbf{k}-\mathbf{p}),[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\perp}(\mathbf{q})]]\\\
&-\frac{1}{3\sqrt{\pi}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{q}\times\mathbf{k}}{k_{\perp}}[b_{\|}(\mathbf{k}-\mathbf{p}),[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\|}(\mathbf{q})]],\\\
\end{split}$ (62)
$\begin{split}\mathfrak{S}_{\eta}^{(5)}(\mathbf{k})=&-\frac{i\mathbf{k}}{\sqrt{\pi}k_{\perp}}\zeta_{(5)}^{i}(\mathbf{k})-i\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}[b_{\|}(\mathbf{k}-\mathbf{p}),-i\mathbf{p}^{i}\zeta_{(3)}^{i}(\mathbf{p})]\\\
&+\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{p}^{l}}{p_{\perp}^{2}}[\zeta_{(3)}^{l}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{p}\times\mathbf{q}}{p_{\perp}^{2}q_{\perp}^{2}}[[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\perp}(\mathbf{q})],b_{\eta}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{p}\cdot\mathbf{q}}{2p_{\perp}^{2}}[[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\|}(\mathbf{q})],b_{\eta}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{2\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}[b_{\|}(\mathbf{k}-\mathbf{p}),[b_{\|}(\mathbf{p}-\mathbf{q}),b_{\eta}(\mathbf{q})]].\\\
\end{split}$ (63)

###  4.2 The bulk terms

In this section we evaluate the bulk terms $\mathfrak{B}_{\eta,\perp}^{(3)}$
and $\mathfrak{B}_{\eta,\perp}^{(5)}$.

#### The bulk terms $\mathfrak{B}_{\eta}^{(3)}$ and
$\mathfrak{B}_{\perp}^{(3)}$

The order-$g^{3}$ bulk terms are defined by

$\begin{split}&\mathfrak{B}_{\eta}^{(3)}(\mathbf{k})=-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{1}^{(2)}(k_{\perp}\tau)S_{\eta}^{(3)}(\tau,\mathbf{k}),\\\
&\mathfrak{B}_{\perp}^{(3)}(\mathbf{k})=-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{0}^{(2)}(k_{\perp}\tau)S_{\perp}^{(3)}(\tau,\mathbf{k}).\\\
\end{split}$ (64)

Here the corresponding source terms are

$\begin{split}S_{\eta}^{(3)}(\tau,\mathbf{k})=&-i\tau^{3}\int
d^{2}\mathbf{x}e^{i\mathbf{k}\cdot\mathbf{x}}\left([\partial_{i}\beta_{i}^{(1)},\beta^{(1)}]+2[\beta_{i}^{(1)},\partial_{i}\beta^{(1)}]\right)\\\
=&i\tau^{2}\int\frac{d^{2}\mathbf{p}}{(2\pi)}\Bigg{(}\frac{2\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|}J_{0}(p_{\perp}\tau)J_{1}(|\mathbf{k}-\mathbf{p}|\tau)\Big{[}b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\Big{]}\Bigg{)}\\\
\end{split}$ (65)

and

$\begin{split}&S_{\perp}^{(3)}(\tau,\mathbf{k})=\frac{-i\epsilon^{ji}\mathbf{k}_{j}}{k_{\perp}}S_{i}^{(3)}\\\
=&-i\tau^{2}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Bigg{(}\frac{\mathbf{p}\times\mathbf{k}}{k_{\perp}p_{\perp}|\mathbf{k}-\mathbf{p}|}\Big{[}b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\Big{]}J_{1}(p_{\perp}\tau)J_{1}(|\mathbf{k}-\mathbf{p}|\tau)\\\
&+\frac{(\mathbf{k}\times\mathbf{p})(\mathbf{k}\cdot\mathbf{p}-p^{2}_{\perp}-k_{\perp}^{2})}{k_{\perp}p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\Big{[}b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\Big{]}J_{0}(p_{\perp}\tau)J_{0}(|\mathbf{k}-\mathbf{p}|\tau)\Bigg{)}\,.\\\
\end{split}$ (66)

From Eq. (65) and Eq. (66) we obtain the order-$g^{3}$ bulk terms

$\begin{split}\mathfrak{B}_{\eta}^{(3)}(\mathbf{k})=&\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left(\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}+i\right)\left[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\right],\\\
\mathfrak{B}_{\perp}^{(3)}(\mathbf{k})=&-\frac{1}{2\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{p}\times\mathbf{k}}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left(\frac{\mathbf{p}\cdot(\mathbf{p}-\mathbf{k})}{|\mathbf{p}\times\mathbf{k}|}-i\right)\left[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})\right]\\\
&+\frac{(\mathbf{p}\cdot(\mathbf{k}-\mathbf{p})-k_{\perp}^{2})}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{k}\times\mathbf{p}}{|\mathbf{p}\times\mathbf{k}|}\left[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})\right]\,.\\\
\end{split}$ (67)

In computing the bulk terms, we need to calculate the proper time integral
from $\tau=0$ to $\tau=\infty$. The integrands involve products of Bessel
functions. The integral formulas used are summarized in Appendix A. Results in
eq. (67) have already been obtained in ref. McLerran:2016snu . We now move on
to compute the new order-$g^{5}$ bulk terms.

#### The bulk term $\mathfrak{B}_{\eta}^{(5)}(\mathbf{k})$

The first task is to compute explicitly the source term at order-$g^{5}$

$\begin{split}S^{(5)}_{\eta}(\tau,\mathbf{k})=&-g\tau^{3}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}(\mathbf{k}+\mathbf{q})_{i}\left[\beta_{i}^{(1)}(\tau,\mathbf{k}-\mathbf{q}),\beta^{(3)}(\tau,\mathbf{q})\right]\\\
&-g\tau^{3}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}(2\mathbf{k}-\mathbf{q})_{i}\left[\beta_{i}^{(3)}(\tau,\mathbf{q}),\beta^{(1)}(\tau,\mathbf{k}-\mathbf{q})\right]\\\
&-g^{2}\tau^{3}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\left[\beta_{i}^{(1)}(\tau,\mathbf{k}-\mathbf{q}),\left[\beta_{i}^{(1)}(\tau,\mathbf{p}),\beta^{(1)}(\tau,\mathbf{q}-\mathbf{p})\right]\right],\end{split}$
(68)

We substitute the solutions $\beta^{(1)}(\tau,\mathbf{k})$,
$\beta_{i}^{(1)}(\tau,\mathbf{k})$, $\beta^{(3)}(\tau,\mathbf{k})$,
$\beta_{i}^{(3)}(\tau,\mathbf{k})$ given in Eqs. (16), (17), (18), and (19) in
order to obtain the explicit expression for $S_{\eta}^{(5)}(\tau,\mathbf{k})$,

$\begin{split}S^{(5)}_{\eta}(\tau,\mathbf{k})=&-i\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{2\mathbf{q}\times\mathbf{k}}{q_{\perp}|\mathbf{k}-\mathbf{q}|^{2}}\left[2\beta^{(3)}(\tau=0,\mathbf{q}),b_{\perp}(\mathbf{k}-\mathbf{q})\right]J_{1}(q_{\perp}\tau)J_{0}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{2\mathbf{q}\times\mathbf{k}}{|\mathbf{k}-\mathbf{q}|^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\Big{[}[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}\\\
&\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})+w_{\perp}^{2}}{q_{\perp}^{2}-w_{\perp}^{2}}\left(\frac{J_{1}(w_{\perp}\tau)}{w_{\perp}}-\frac{J_{1}(q_{\perp}\tau)}{q_{\perp}}\right)J_{0}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+i\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{2\mathbf{k}\times\mathbf{q}}{q_{\perp}|\mathbf{k}-\mathbf{q}|}\Big{[}\beta_{\perp}^{(3)}(\tau=0,\mathbf{q}),b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}J_{0}(q_{\perp}\tau)J_{1}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{2\mathbf{k}\times\mathbf{q}}{q^{2}_{\perp}|\mathbf{k}-\mathbf{q}|}\frac{\mathbf{q}\times\mathbf{p}}{2p^{2}_{\perp}|\mathbf{q}-\mathbf{p}|^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{q_{\perp}^{2}-w_{\perp}^{2}}\Big{(}(-p_{\perp}^{2}-|\mathbf{q}-\mathbf{p}|^{2}+w_{\perp}^{2})\\\
&\qquad\Big{[}\big{[}b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\big{]},b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}+2(\mathbf{p}\cdot\mathbf{q}-p_{\perp}^{2}-q_{\perp}^{2})\Big{[}\big{[}b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})\big{]},b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\Big{)}\\\
&\qquad\qquad\times(J_{0}(w_{\perp}\tau)-J_{0}(q_{\perp}\tau))J_{1}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(2\mathbf{k}-\mathbf{q})\cdot\mathbf{q}}{q_{\perp}^{2}|\mathbf{k}-\mathbf{q}|}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})}{4p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}}\Big{(}(-q_{\perp}^{2}+w_{\perp}^{2}+2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p}))\\\
&\qquad\Big{[}\left[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\right],b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}+2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\Big{[}\left[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})\right],b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\Big{)}\\\
&\qquad\qquad\times(1-J_{0}(w_{\perp}\tau))J_{1}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(\mathbf{k}-\mathbf{q})\cdot\mathbf{p}}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}||\mathbf{k}-\mathbf{q}|^{2}}\Big{[}b_{\perp}(\mathbf{k}-\mathbf{q}),[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})]\Big{]}\\\
&\qquad\qquad\times
J_{1}(|\mathbf{q}-\mathbf{p}|\tau)J_{0}(p_{\perp}\tau)J_{0}(|\mathbf{k}-\mathbf{q}|\tau).\\\
\end{split}$ (69)

Using the above result, we can evaluate the bulk term
$\mathfrak{B}_{\eta}^{(5)}$ by

$\begin{split}\mathfrak{B}_{\eta}^{(5)}(\mathbf{k})=-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{1}^{(2)}(k_{\perp}\tau)S_{\eta}^{(5)}(\tau,\mathbf{k}).\end{split}$
(70)

Again, it requires integrations over the proper time from $\tau=0$ to
$\tau=\infty$. 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

$\begin{split}\mathfrak{B}^{(5)}_{\eta}(\mathbf{k})=&\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\left[2\beta^{(3)}(\tau=0,\mathbf{q}),b_{\perp}(\mathbf{k}-\mathbf{q})\right]\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}|\mathbf{k}-\mathbf{q}|^{2}}\left(\frac{\mathbf{k}\cdot\mathbf{q}}{|\mathbf{k}\times\mathbf{q}|}+i\right)\\\
&+\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\Big{[}q_{\perp}\beta_{\perp}^{(3)}(\tau=0,\mathbf{q}),b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}|\mathbf{k}-\mathbf{q}|^{2}}\left(\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}+i\right)\\\
&+\frac{i}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{[}[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathcal{I}_{1}(\mathbf{p},\mathbf{q},\mathbf{k})}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}}\\\
&+\frac{i}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{[}\big{[}b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\big{]},b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathcal{I}_{2}(\mathbf{p},\mathbf{q},\mathbf{k})}{p^{2}_{\perp}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}}\\\
&+\frac{i}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{[}\big{[}b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})\big{]},b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathcal{I}_{3}(\mathbf{p},\mathbf{q},\mathbf{k})}{p^{2}_{\perp}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}}\,.\\\
\end{split}$ (71)

Here $\mathcal{I}_{1,2,3}(\mathbf{p},\mathbf{q},\mathbf{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:

$\begin{split}&\mathcal{I}_{1}(\mathbf{p},\mathbf{q},\mathbf{k})\\\ =&\pi
k_{\perp}(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{|\mathbf{q}-\mathbf{p}|^{2}-p_{\perp}^{2}+w_{\perp}^{2}}{q_{\perp}^{2}-w_{\perp}^{2}}\\\
&\qquad\times\Big{(}\frac{1}{w_{\perp}}L_{011}(|\mathbf{k}-\mathbf{q}|,w_{\perp},k_{\perp})-\frac{1}{q_{\perp}}L_{011}(|\mathbf{k}-\mathbf{q}|,q_{\perp},k_{\perp})\Big{)}\\\
&+\frac{1}{4}\pi
k_{\perp}(\mathbf{k}-\mathbf{q})\cdot\mathbf{p}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{|\mathbf{q}-\mathbf{p}|^{2}-p_{\perp}^{2}+w_{\perp}^{2}}{w_{\perp}}L_{011}(|\mathbf{k}-\mathbf{q}|,w_{\perp},k_{\perp}),\\\
\end{split}$ (72)
$\begin{split}&\mathcal{I}_{2}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
=&\frac{1}{2q_{\perp}^{2}}(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})\pi|\mathbf{k}-\mathbf{q}|k_{\perp}\\\
&\quad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\left(1+\frac{2\mathbf{p}\cdot(\mathbf{p}-\mathbf{q})}{q_{\perp}^{2}-w_{\perp}^{2}}\right)\Big{(}L_{011}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-L_{011}(q_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)}\\\
&+\frac{1}{8q_{\perp}^{2}}(k_{\perp}^{2}-|\mathbf{k}-\mathbf{q}|^{2})(|\mathbf{q}-\mathbf{p}|^{2}-p_{\perp}^{2})\pi|\mathbf{k}-\mathbf{q}|k_{\perp}\\\
&\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\left(1-\frac{|\mathbf{q}-\mathbf{p}|^{2}+p_{\perp}^{2}}{w_{\perp}^{2}}\right)\Big{(}L_{011}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-L_{011}(0,|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)},\\\
\end{split}$ (73)
$\begin{split}&\mathcal{I}_{3}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
=&\frac{1}{q_{\perp}^{2}}(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})(p_{\perp}^{2}+q_{\perp}^{2}-\mathbf{p}\cdot\mathbf{q})\pi|\mathbf{k}-\mathbf{q}|k_{\perp}\\\
&\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{q_{\perp}^{2}-w_{\perp}^{2}}\Big{(}L_{011}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-L_{011}(q_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)}\\\
&+\frac{1}{4q_{\perp}^{2}}\left(k_{\perp}^{2}-|\mathbf{k}-\mathbf{q}|^{2}\right)\left(|\mathbf{q}-\mathbf{p}|^{2}-p_{\perp}^{2}\right)\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\pi|\mathbf{k}-\mathbf{q}|k_{\perp}\\\
&\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}}\Big{(}L_{011}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-L_{011}(0,|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)},\\\
\end{split}$ (74)

where
$w_{\perp}=\sqrt{p_{\perp}^{2}+|\mathbf{q}-\mathbf{p}|^{2}-2p_{\perp}|\mathbf{q}-\mathbf{p}|\cos\phi}$.
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
$\mathcal{I}_{1,2,3}(\mathbf{p},\mathbf{q},\mathbf{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 $\mathfrak{B}_{\perp}^{(5)}$

For the transverse polarization, the corresponding order-$g^{5}$ source term
has the expression

$\begin{split}S^{(5)}_{\perp}(\tau,\mathbf{k})=&\frac{i\epsilon^{il}\mathbf{k}_{l}}{k_{\perp}}\tau^{2}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}(\mathbf{k}-2\mathbf{q})_{i}\left([\tilde{\beta}^{(3)}(\tau,\mathbf{q}),\tilde{\beta}^{(1)}(\tau,\mathbf{k}-\mathbf{q})]+[\beta_{j}^{(3)}(\tau,\mathbf{q}),\beta_{j}^{(1)}(\tau,\mathbf{k}-\mathbf{q})]\right)\\\
&\qquad+\left(-(2\mathbf{k}-\mathbf{q})_{j}[\beta_{j}^{(3)}(\tau,\mathbf{q}),\beta_{i}^{(1)}(\tau,\mathbf{k}-\mathbf{q})]+(\mathbf{k}+\mathbf{q})_{j}[\beta_{i}^{(3)}(\tau,\mathbf{q}),\beta_{j}^{(1)}(\tau,\mathbf{k}-\mathbf{q})]\right)\\\
&\qquad+\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{(}[\tilde{\beta}^{(1)}(\tau,\mathbf{k}-\mathbf{q}),[\beta_{i}^{(1)}(\tau,\mathbf{p}),\tilde{\beta}^{(1)}(\tau,\mathbf{q}-\mathbf{p})]]\\\
&\qquad\qquad+[\beta_{j}^{(1)}(\tau,\mathbf{k}-\mathbf{q}),[\beta_{i}^{(1)}(\tau,\mathbf{p}),\beta_{j}^{(1)}(\tau,\mathbf{q}-\mathbf{p})]]\Big{)}\,.\end{split}$
(75)

Again, using the solutions $\beta^{(1)}(\tau,\mathbf{k})$,
$\beta_{i}^{(1)}(\tau,\mathbf{k})$, $\beta^{(3)}(\tau,\mathbf{k})$,
$\beta_{i}^{(3)}(\tau,\mathbf{k})$ given in Eqs. (16), (17) (18), and (19),
one obtains its explicit expression

$\begin{split}&S_{\perp}^{(5)}(\tau,\mathbf{k})\\\
=&i\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{2(\mathbf{k}\times\mathbf{q})}{q_{\perp}|\mathbf{k}-\mathbf{q}|}\Big{[}2\beta^{(3)}(\tau=0,\mathbf{q}),b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}J_{1}(q_{\perp}\tau)J_{1}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&-\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{2(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|}\Big{[}[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\\\
&\qquad\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})+w_{\perp}^{2}}{q_{\perp}^{2}-w_{\perp}^{2}}\left(\frac{J_{1}(w_{\perp}\tau)}{w_{\perp}}-\frac{J_{1}(q_{\perp}\tau)}{q_{\perp}}\right)J_{1}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+i\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})(2\mathbf{k}\times\mathbf{q})}{q_{\perp}|\mathbf{k}-\mathbf{q}|^{2}}[\beta_{\perp}^{(3)}(\tau=0,\mathbf{q}),b_{\perp}(\mathbf{k}-\mathbf{q})]J_{0}(q_{\perp}\tau)J_{0}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})}{q^{2}_{\perp}|\mathbf{k}-\mathbf{q}|^{2}p^{2}_{\perp}|\mathbf{q}-\mathbf{p}|^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{q_{\perp}^{2}-w_{\perp}^{2}}\Big{(}(w_{\perp}^{2}-p_{\perp}^{2}-|\mathbf{q}-\mathbf{p}|^{2})\\\
&\qquad\qquad\times\Big{[}[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}+2(\mathbf{p}\cdot\mathbf{q}-p_{\perp}^{2}-q_{\perp}^{2})\Big{[}[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}\Big{)}\\\
&\qquad\qquad\times(J_{0}(w_{\perp}\tau)-J_{0}(q_{\perp}\tau))J_{0}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(2\mathbf{k}-\mathbf{q})\cdot\mathbf{q}\,\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})}{4p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}q_{\perp}^{2}}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}}\Big{(}(w_{\perp}^{2}-q_{\perp}^{2}+2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p}))\\\
&\qquad\qquad\times\left[[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\right]+2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\Big{[}b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}\Big{)}\\\
&\qquad\qquad\times(1-J_{0}(w_{\perp}\tau))J_{0}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&+\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{p}\cdot\mathbf{k}}{|\mathbf{q}-\mathbf{p}||\mathbf{k}-\mathbf{q}|p_{\perp}^{2}}\Big{[}b_{\eta}(\mathbf{k}-\mathbf{q}),[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})]\Big{]}\\\
&\qquad\quad\times
J_{1}(|\mathbf{q}-\mathbf{p}|\tau)J_{0}(p_{\perp}\tau)J_{1}(|\mathbf{k}-\mathbf{q}|\tau)\\\
&-\tau^{2}\frac{1}{k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(\mathbf{k}-\mathbf{q})\cdot(\mathbf{q}-\mathbf{p})\mathbf{p}\cdot\mathbf{k}}{|\mathbf{q}-\mathbf{p}|^{2}p_{\perp}^{2}|\mathbf{k}-\mathbf{q}|^{2}}\Big{[}b_{\perp}(\mathbf{k}-\mathbf{q}),[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})]\Big{]}\\\
&\qquad\quad\times
J_{0}(|\mathbf{q}-\mathbf{p}|\tau)J_{0}(p_{\perp}\tau)J_{0}(|\mathbf{k}-\mathbf{q}|\tau).\\\
\end{split}$ (76)

The source term is then used to compute the order-$g^{5}$ bulk term,
contributing to single gluon production amplitude, by

$\mathfrak{B}_{\perp}^{(5)}(\mathbf{p})=-i\sqrt{\frac{\pi}{4}}\int_{0}^{\infty}d\tau\frac{1}{\tau}H_{0}^{(2)}(p_{\perp}\tau)S^{(5)}_{\perp}(\tau,\mathbf{p}).$
(77)

We skip the details of doing the proper time integral and present the final
result here

$\begin{split}\mathfrak{B}_{\perp}^{(5)}(\mathbf{k})=&\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\Big{[}2\beta^{(3)}(\tau=0,\mathbf{q}),b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}|\mathbf{k}-\mathbf{q}|^{2}}\left(\frac{\mathbf{q}\cdot(\mathbf{q}-\mathbf{k})}{|\mathbf{q}\times\mathbf{k}|}-i\right)\\\
&+\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\left[q_{\perp}\beta_{\perp}^{(3)}(\tau=0,\mathbf{q}),b_{\perp}(\mathbf{k}-\mathbf{q})\right]\frac{(\mathbf{k}\times\mathbf{q})(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})}{q^{2}_{\perp}|\mathbf{k}-\mathbf{q}|^{2}|\mathbf{k}\times\mathbf{q}|}\\\
&+\frac{i}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{[}\left[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})\right],b_{\eta}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathcal{G}_{1}(\mathbf{p},\mathbf{q},\mathbf{k})}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}}\\\
&+i\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{[}[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathcal{G}_{2}(\mathbf{p},\mathbf{q},\mathbf{k})}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}}\\\
&+i\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\Big{[}[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{q}-\mathbf{p})],b_{\perp}(\mathbf{k}-\mathbf{q})\Big{]}\frac{\mathcal{G}_{3}(\mathbf{p},\mathbf{q},\mathbf{k})}{p_{\perp}^{2}|\mathbf{q}-\mathbf{p}|^{2}|\mathbf{k}-\mathbf{q}|^{2}}.\\\
\end{split}$ (78)

Similar to the analysis for $\mathfrak{B}_{\eta}^{(5)}(\mathbf{k})$ in the
previous section, we introduced another three auxilliary functions
$\mathcal{G}_{1,2,3}(\mathbf{p},\mathbf{q},\mathbf{k})$ defined as

$\begin{split}&\mathcal{G}_{1}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
=&\pi|\mathbf{k}-\mathbf{q}|(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})+w_{\perp}^{2}}{q_{\perp}^{2}-w_{\perp}^{2}}\left(\frac{1}{w_{\perp}}L_{110}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-\frac{1}{q_{\perp}}L_{110}(q_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})\right)\\\
&+\pi|\mathbf{k}-\mathbf{q}|(\mathbf{p}\cdot\mathbf{k})\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{|\mathbf{q}-\mathbf{p}|^{2}-p_{\perp}^{2}+w_{\perp}^{2}}{4w_{\perp}}L_{110}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp}),\\\
\end{split}$ (79)
$\begin{split}\mathcal{G}_{2}(\mathbf{p},\mathbf{q},\mathbf{k})=&\frac{1}{q_{\perp}^{2}}(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})\frac{\pi}{2}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{q_{\perp}^{2}-w_{\perp}^{2}}(-w_{\perp}^{2}+p_{\perp}^{2}+|\mathbf{q}-\mathbf{p}|^{2})\\\
&\times\Big{(}L_{000}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-L_{000}(q_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)}\\\
&+\frac{1}{4q_{\perp}^{2}}(2\mathbf{k}-\mathbf{q})\cdot\mathbf{q}\,\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})\frac{\pi}{2}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}}(w_{\perp}^{2}-q_{\perp}^{2}+2\mathbf{p}\cdot(\mathbf{q}-\mathbf{p}))\\\
&\times\Big{(}L_{000}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp}-L_{000}(0,|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)}\,,\\\
\end{split}$ (80)
$\begin{split}\mathcal{G}_{3}(\mathbf{p},\mathbf{q},\mathbf{k})=&\frac{\pi}{q_{\perp}^{2}}(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})(-\mathbf{p}\cdot\mathbf{q}+p_{\perp}^{2}+q_{\perp}^{2})\\\
&\qquad\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{q_{\perp}^{2}-w_{\perp}^{2}}\Big{(}L_{000}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})-L_{000}(q_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)}\\\
&+\frac{\pi}{4q_{\perp}^{2}}(2\mathbf{k}-\mathbf{q})\cdot\mathbf{q}\,\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})\mathbf{q}\cdot(\mathbf{q}-2\mathbf{p})\mathbf{p}\cdot(\mathbf{q}-\mathbf{p})\\\
&\qquad\qquad\times\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}\frac{1}{w_{\perp}^{2}}\Big{(}L_{000}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp}-L_{000}(0,|\mathbf{k}-\mathbf{q}|,k_{\perp})\Big{)}\\\
&-\frac{1}{2}(\mathbf{k}\times\mathbf{q})(\mathbf{q}\times\mathbf{p})\frac{\pi}{2}\int_{-\pi}^{\pi}\frac{d\phi}{2\pi}L_{000}(w_{\perp},|\mathbf{k}-\mathbf{q}|,k_{\perp}).\\\
\end{split}$ (81)

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
$\mathbf{p},\mathbf{q},\mathbf{k}$.

##  5 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 $\rho_{P}^{a}(\mathbf{k})$ and the target Wilson
line $U^{ab}(\mathbf{k})$. The following explicit expressions will be
repeatedly used. We list them for reference.

$\begin{split}&b^{c}_{\eta}(\mathbf{k})=\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}}{p_{\perp}^{2}}\rho_{P}^{a}(\mathbf{p})U^{ac}(\mathbf{k}-\mathbf{p}),\\\
&b^{c}_{\perp}(\mathbf{k})=\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\rho_{P}^{a}(\mathbf{p})U^{ac}(\mathbf{k}-\mathbf{p})\,,\\\
&b_{\|}^{c}(\mathbf{k})=\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{-\mathbf{k}\cdot\mathbf{p}}{k_{\perp}^{2}p_{\perp}^{2}}\rho_{P}^{a}(\mathbf{p})U^{ac}(\mathbf{k}-\mathbf{p})\,.\\\
\end{split}$ (82)

Further more, the eikonally color rotated order-$g^{3}$ and order-$g^{5}$
Weiszacker-Williams fields are

$\begin{split}&\zeta_{(3)}^{c,l}(\mathbf{q})=\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\alpha_{P,(3)}^{a,l}(\mathbf{p})U^{ac}(\mathbf{q}-\mathbf{p}),\\\
&\zeta_{(5)}^{c,l}(\mathbf{k})=\int\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\alpha_{P,(5)}^{a,l}(\mathbf{q})U^{ac}(\mathbf{k}-\mathbf{q})\\\
\end{split}$ (83)

with

$\begin{split}\alpha_{P,(3)}^{a,l}(\mathbf{p})=&\frac{i}{2}\left(\delta^{lj}-\frac{\mathbf{p}^{l}\mathbf{p}^{j}}{p_{\perp}^{2}}\right)\int\frac{d^{2}\mathbf{p}_{1}}{(2\pi)^{2}}\frac{-1}{|\mathbf{p}-\mathbf{p}_{1}|^{2}}\frac{i\mathbf{p}_{1}^{j}}{p_{1}^{2}}\Big{[}\rho_{P}(\mathbf{p}-\mathbf{p}_{1}),\rho_{P}(\mathbf{p}_{1})\Big{]}\\\
=&\frac{1}{2}\left(\delta^{lj}-\frac{\mathbf{p}^{l}\mathbf{p}^{j}}{p_{\perp}^{2}}\right)\int\frac{d^{2}\mathbf{p}_{1}}{(2\pi)^{2}}\frac{\mathbf{p}_{1}^{j}}{|\mathbf{p}-\mathbf{p}_{1}|^{2}p_{1}^{2}}\Big{[}\rho_{P}(\mathbf{p}-\mathbf{p}_{1}),\rho_{P}(\mathbf{p}_{1})\Big{]},\\\
\end{split}$ (84)
$\begin{split}\alpha_{P,(5)}^{a,l}(\mathbf{q})=&\left(\delta^{ij}-\frac{\mathbf{q}^{i}\mathbf{q}^{j}}{q_{\perp}^{2}}\right)\Big{(}\frac{1}{2}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}_{1}}{(2\pi)^{2}}\frac{1}{p_{\perp}^{2}p_{1}^{2}}\frac{i(\mathbf{q}-\mathbf{p})^{j}}{|\mathbf{q}-\mathbf{p}|^{2}}[[\rho_{P}(\mathbf{p}_{1}),\rho_{P}(\mathbf{p}-\mathbf{p}_{1})],\rho_{P}(\mathbf{q}-\mathbf{p})]\\\
&+\frac{1}{6}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{p}_{1}}{(2\pi)^{2}}\frac{i(\mathbf{p}-\mathbf{p}_{1})^{j}}{|\mathbf{p}-\mathbf{p}_{1}|^{2}}\frac{1}{p_{1}^{2}|\mathbf{q}-\mathbf{p}|^{2}}[[\rho_{P}(\mathbf{p}_{1}),\rho_{P}(\mathbf{p}-\mathbf{p}_{1})],\rho_{P}(\mathbf{q}-\mathbf{p})]\Big{)}.\end{split}$
(85)

###  5.1 Leading order gluon production

The leading order gluon production is defined by

$n^{(2)}(\mathbf{k})=|\mathfrak{S}_{\eta}^{(1)}(\mathbf{k})|^{2}+|\mathfrak{S}_{\perp}^{(1)}(\mathbf{k})|^{2}\,.$
(86)

Using the explicit forms for

$\mathfrak{S}_{\eta}^{(1)}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}b_{\eta}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}}{p_{\perp}^{2}}\rho^{a}(\mathbf{p})U^{ac}(\mathbf{k}-\mathbf{p})T^{c},$
(87)
$\mathfrak{S}_{\perp}^{(1)}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}b_{\perp}(\mathbf{k})=\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\rho^{a}(\mathbf{p})U^{ac}(\mathbf{k}-\mathbf{p})T^{c}$
(88)

we obtain that at the leading order

$\begin{split}&\frac{dN_{\rm
LO}}{d^{2}\mathbf{k}dy}=\frac{1}{(2\pi)^{2}}g^{2}n^{(2)}(\mathbf{k})\\\
=&\frac{g^{2}}{(2\pi)^{3}}\frac{1}{\pi
k_{\perp}^{2}}(\delta^{ij}\delta^{mn}+\epsilon^{ij}\epsilon^{mn})\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}\frac{\mathbf{p}_{i}(\mathbf{k}-\mathbf{p})_{j}}{\mathbf{p}^{2}}\frac{\mathbf{q}_{m}(\mathbf{k}-\mathbf{q})_{n}}{\mathbf{q}^{2}}\\\
&\times\rho^{a}(\mathbf{p})\left[U(\mathbf{k}-\mathbf{p})U^{\dagger}(\mathbf{k}-\mathbf{q})\right]_{ad}\rho^{d}(-\mathbf{q})\,.\end{split}$
(89)

This result reproduces the dilute-dense single gluon spectrum calculated in
Refs. Dumitru:2001ux ; Blaizot:2004wu ; Blaizot:2008yb . Note that this
expression is given in the functional form with explicit dependence on
$\rho_{P}$ and $\rho_{T}$.

For comparison with high order corrections, we use the Lipatov vertex
Fadin:1975cb ; Kuraev:1977fs ; Kovchegov:2012mbw

$L^{i}(\mathbf{p},\mathbf{k})=\left(\frac{\mathbf{p}^{i}}{p_{\perp}^{2}}-\frac{\mathbf{k}^{i}}{k_{\perp}^{2}}\right).$
(90)

Then the leading order gluon production can be expressed in a more informative
form

$n^{(2)}(\mathbf{k})=-\frac{1}{\pi}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{d^{2}\mathbf{q}}{(2\pi)^{2}}L^{i}(\mathbf{p},\mathbf{k})L^{i}(\mathbf{q},-\mathbf{k})\rho^{a}(\mathbf{p})\rho^{b}(\mathbf{q})U^{ac}(\mathbf{k}-\mathbf{p})U^{bc}(-\mathbf{k}-\mathbf{q}).$
(91)

The color structure can be shown by Fig. 4. The coefficient function is just
the Lipatov vertex squared.

Figure 4: Color structure for $|M_{(1)}|^{2}$. Two gluons lines cross the
target nucleus indicating the two Wilson lines. The dashed line separates the
amplitude on the left and the complex conjugate amplitude on the right.

###  5.2 Next-to-leading order gluon production

The order-$g^{4}$ contribution to the single gluon production is defined by

$n^{(4)}(\mathbf{k})=\sum_{\gamma=\eta,\perp}\Big{[}\mathfrak{S}_{\gamma}^{(1)\ast}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(3)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)}(\mathbf{k})\right)+\mathfrak{S}_{\gamma}^{(1)}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(3)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)\ast}(\mathbf{k})\right)\Big{]}\,.$
(92)

The surface terms $\mathfrak{S}_{\perp}^{(1)}(\mathbf{k})$,
$\mathfrak{S}_{\eta}^{(1)}(\mathbf{k})$,
$\mathfrak{S}_{\perp}^{(3)}(\mathbf{k})$ and
$\mathfrak{S}_{\eta}^{(3)}(\mathbf{k})$ have been calculated in Eqs. (59),
(60), and (61). The bulk terms at order-$g^{3}$ are given in Eqs. (67).

Collecting the results for the surface and bulk terms at order-$g^{3}$, we
have

$\begin{split}&\mathfrak{S}_{\eta}^{(3)}(\mathbf{k})+\mathfrak{B}_{\eta}^{(3)}(\mathbf{k})\\\
=&\frac{-i\mathbf{k}^{i}}{\sqrt{\pi}k_{\perp}}\zeta_{(3)}^{i}(\mathbf{k})+\frac{ig}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)}\left[b_{\eta}(\mathbf{p}),b_{\|}(\mathbf{k}-\mathbf{p})\right]\\\
&+\frac{1}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left(\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}+i\right)[b_{\perp}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})],\\\
\end{split}$ (93)

and

$\begin{split}&\mathfrak{S}_{\perp}^{(3)}(\mathbf{k})+\mathfrak{B}_{\perp}^{(3)}(\mathbf{k})\\\
=&\frac{i\epsilon^{ih}\mathbf{k}^{h}}{\sqrt{\pi}k_{\perp}}\zeta_{(3)}^{i}(\mathbf{k})+\frac{ig}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{k}\cdot\mathbf{p}}{p_{\perp}^{2}}[b_{\perp}(\mathbf{p}),b_{\|}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{2}\frac{ig}{\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}(\mathbf{k}\times\mathbf{p})[b_{\|}(\mathbf{p}),b_{\|}(\mathbf{k}-\mathbf{p})]\\\
&-\frac{1}{2\sqrt{\pi}k_{\perp}}\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}\frac{\mathbf{p}\times\mathbf{k}}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left(\frac{\mathbf{p}\cdot(\mathbf{p}-\mathbf{k})}{|\mathbf{p}\times\mathbf{k}|}-i\right)[b_{\eta}(\mathbf{p}),b_{\eta}(\mathbf{k}-\mathbf{p})]\\\
&+\frac{(\mathbf{p}\cdot(\mathbf{k}-\mathbf{p})-k_{\perp}^{2})}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{k}\times\mathbf{p}}{|\mathbf{p}\times\mathbf{k}|}[b_{\perp}(\mathbf{p}),b_{\perp}(\mathbf{k}-\mathbf{p})]\,.\end{split}$
(94)

From these surface terms and bulk terms, after some tediou algebra,
substituting eqs (82), one obtains the next-to-leading order single inclusive
gluon production

$\begin{split}n^{(4)}(\mathbf{k})=&\frac{1}{\pi}\int_{\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2}}\mathcal{E}_{1}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{k})f^{ade}\rho_{P}^{d}(\mathbf{p}-\mathbf{p}_{1})\rho_{P}^{e}(\mathbf{p}_{1})\rho_{P}^{e_{2}}(\mathbf{p}_{2})U^{ac}(\mathbf{k}-\mathbf{p})U^{e_{2}c}(-\mathbf{k}-\mathbf{p}_{2})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3}}\mathcal{E}_{2}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})f^{cde}\rho_{P}^{e_{1}}(\mathbf{p}_{1})\rho_{P}^{e_{2}}(\mathbf{p}_{2})\rho_{P}^{e_{3}}(\mathbf{p}_{3})\\\
&\qquad\times
U^{e_{1}d}(\mathbf{p}-\mathbf{p}_{1})U^{e_{2}e}(\mathbf{k}-\mathbf{p}-\mathbf{p}_{2})U^{e_{3}c}(-\mathbf{k}-\mathbf{p}_{3})\\\
&+c.c.\end{split}$ (95)

Here $c.c.$ represents the complex conjugate terms and we have introduced the
shorthand notation $\int_{\mathbf{p}}=\int\frac{d^{2}\mathbf{p}}{(2\pi)^{2}}$.
The color structure (functional dependence on $\rho_{P}$ and $U$) are
correctly captured by the two diagrams in Fig. 5. The kinematic function
$\mathcal{E}_{1}$ has the expression

$\mathcal{E}_{1}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{k})=\frac{1}{2}\frac{(\mathbf{p}_{1}\times\mathbf{p})}{|\mathbf{p}-\mathbf{p}_{1}|^{2}p_{1}^{2}}\left(\frac{(\mathbf{p}\times\mathbf{p}_{2})}{p_{\perp}^{2}p_{2}^{2}}+\frac{(\mathbf{p}\times\mathbf{k})}{p_{\perp}^{2}k^{2}_{\perp}}\right).$
(96)

One immediately sees that the complex conjugate amplitude is represented by
the Lipatov vertex

$L^{j}(\mathbf{p}_{2},-\mathbf{k})=\left(\frac{\mathbf{p}_{2}^{j}}{p_{2}^{2}}+\frac{\mathbf{k}^{j}}{k_{\perp}^{2}}\right).$
(97)

We thus introduce an effective vertex at order-$g^{3}$,

$\Gamma^{j}_{1}(\mathbf{p},\mathbf{p}_{1})=\frac{1}{2}\frac{(\mathbf{p}_{1}\times\mathbf{p})}{|\mathbf{p}-\mathbf{p}_{1}|^{2}p_{1}^{2}}\frac{\epsilon^{ij}\mathbf{p}^{i}}{p_{\perp}^{2}}.$
(98)

The coefficient function in eq. (96) follows simply as the product of the two
effective vertices.

$\mathcal{E}_{1}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{k})=\Gamma^{j}_{1}(\mathbf{p}_{1},\mathbf{p})L^{j}(\mathbf{p}_{2},-\mathbf{k})$
(99)

See Fig. 5(a) for the momentum assignments and the color structure.

(a)

(b)

Figure 5: Color structure for terms in the $M^{\ast}_{(1)}M_{(3)}$ that
contain two and three Wilson lines: (a) two gluon lines cross the target
nucleus and (b) three gluon lines cross the target nucleus. The dashed line
separates the amplitude on the left and the complex conjugate amplitude on the
right.

The second kinematic function in eq. (95) is

$\begin{split}&\mathcal{E}_{2}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\\\
=&-\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{k}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}+\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{k^{2}_{\perp}p_{3}^{2}}\\\
&-\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{k}\times\mathbf{p})}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left(i\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}-1\right)\frac{(\mathbf{k}+\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{k^{2}_{\perp}p_{3}^{2}}\\\
&-\left[\frac{(\mathbf{k}\cdot\mathbf{p})(\mathbf{p}\times\mathbf{p}_{1})}{p_{\perp}^{2}p_{1}^{2}}+\frac{1}{2}\frac{(\mathbf{k}\times\mathbf{p})(\mathbf{p}\cdot\mathbf{p}_{1})}{p_{\perp}^{2}p_{1}^{2}}\right]\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{k}-\mathbf{p}|^{2}p_{2}^{2}}\frac{\mathbf{k}\times\mathbf{p}_{3}}{k^{2}_{\perp}p_{3}^{2}}\\\
&-\frac{1}{2}\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{k}\times\mathbf{p})}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left(i\frac{\mathbf{p}\cdot(\mathbf{p}-\mathbf{k})}{|\mathbf{p}\times\mathbf{k}|}+1\right)\frac{\mathbf{k}\times\mathbf{p}_{3}}{k^{2}_{\perp}p_{3}^{2}}\\\
&+\frac{1}{2}\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p})\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{p}\cdot(\mathbf{k}-\mathbf{p})-k_{\perp}^{2})}{p_{\perp}^{2}|\mathbf{k}-\mathbf{p}|^{2}}i\frac{\mathbf{k}\times\mathbf{p}}{|\mathbf{k}\times\mathbf{p}|}\frac{\mathbf{k}\times\mathbf{p}_{3}}{k^{2}_{\perp}p_{3}^{2}}\end{split}$
(100)

Again, isolating the Lipatov vertex in the complex conjugate amplitude

$L^{j}(\mathbf{p}_{3},-\mathbf{k})=\left(\frac{\mathbf{p}_{3}^{j}}{p_{3}^{2}}+\frac{\mathbf{k}^{j}}{k_{\perp}^{2}}\right),$
(101)

we introduce another effective vertex at order-$g^{3}$,

$\begin{split}&\Gamma^{j}_{2}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{k})\\\
=&-\Bigg{\\{}\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{k}-\mathbf{p}|^{2}p_{2}^{2}}+\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{|\mathbf{k}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}\times\mathbf{p})}{p_{\perp}^{2}}\left(i\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}-1\right)\Bigg{\\}}\frac{\mathbf{k}^{j}}{k_{\perp}^{2}}\\\
&-\Bigg{\\{}\left[\frac{(\mathbf{p}\times\mathbf{p}_{1})}{p_{1}^{2}}\frac{(\mathbf{k}\cdot\mathbf{p})}{p_{\perp}^{2}}+\frac{1}{2}\frac{(\mathbf{p}\cdot\mathbf{p}_{1})}{p_{1}^{2}}\frac{(\mathbf{k}\times\mathbf{p})}{p_{\perp}^{2}}\right]\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{k}-\mathbf{p}|^{2}p_{2}^{2}}\\\
&\qquad+\frac{1}{2}\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{(\mathbf{k}\times\mathbf{p})}{p_{\perp}^{2}}\left(i\frac{\mathbf{p}\cdot(\mathbf{p}-\mathbf{k})}{|\mathbf{p}\times\mathbf{k}|}+1\right)\\\
&\qquad+\frac{i}{2}\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{k}-\mathbf{p})\times\mathbf{p}_{2}}{p_{2}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\frac{k_{\perp}^{2}+p_{\perp}^{2}-\mathbf{p}\cdot\mathbf{k}}{|\mathbf{k}\times\mathbf{p}|}\Bigg{\\}}\frac{\epsilon^{lj}\mathbf{k}^{l}}{k_{\perp}^{2}}\,.\end{split}$
(102)

As a consequence, the coefficient function $\mathcal{E}_{2}$ can be expressed
as a product of the two effective vertices

$\mathcal{E}_{2}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})=\Gamma^{j}_{2}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{k})L^{j}(\mathbf{p}_{3},-\mathbf{k})\,.$
(103)

See Fig. 5(b) for the assignments of the various momenta.

Note that $\mathcal{E}_{2}$ has both real and imaginary parts. The imaginary
part is always accompanied by the sign function
$\mathrm{sgn}(\mathbf{k}\times\mathbf{p})=(\mathbf{k}\times\mathbf{p})/|\mathbf{k}\times\mathbf{p}|$.
It is precisely this imaginary part that will be responsible for the asymmetry
between $n^{(4)}(\mathbf{k})$ and $n^{(4)}(-\mathbf{k})$. In
$n^{(4)}(\mathbf{k})$, the odd part under
$\mathbf{k}\leftrightarrow-\mathbf{k}$ has already been extracted in Ref.
McLerran:2016snu . The results were contrasted with those obtained from a
diagrammatic approach of Ref. Chirilli:2015tea  in Ref. Kovchegov:2018jun .
The even part under $\mathbf{k}\leftrightarrow-\mathbf{k}$ is computed here
for the first time. Also note that since $n^{(4)}(\mathbf{k})$ is cubic in
$\rho_{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)}(\mathbf{k})n^{(4)}(\mathbf{p})$ is proportional to $\rho_{P}^{6}$,
which does not vanish in the McLerran-Venugopalan model).

Considering the color structures for $n^{(4)}(\mathbf{k})$ in Eq. (95), 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}(\mathbf{p}_{3},\mathbf{k})$,
$\Gamma_{1}^{j}(\mathbf{p},\mathbf{p}_{1})$,
$\Gamma_{2}^{j}(\mathbf{p},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{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(\mathbf{x})=1$, meaning the absence of the nucleus, there will be no gluon
productions. Using
$U^{ab}(\mathbf{p})\rightarrow\delta^{ab}\delta^{(2)}(\mathbf{p})$, it is easy
to check that $n^{(4)}(\mathbf{k})$ expressed in Eq. (95) vanishes.

###  5.3 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 $\rho_{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

$\begin{split}n^{(6)}(\mathbf{k})=&\sum_{\gamma=\eta,\perp}\Big{[}(\mathfrak{S}_{\gamma}^{(3)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)}(\mathbf{k}))(\mathfrak{S}_{\gamma}^{(3)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)\ast}(\mathbf{k}))\\\
&+\mathfrak{S}_{\gamma}^{(1)\ast}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(5)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(5)}(\mathbf{k})\right)+\mathfrak{S}_{\gamma}^{(1)}(\mathbf{k})\left(\mathfrak{S}_{\gamma}^{(5)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(5)\ast}(\mathbf{k})\right)\Big{]}.\\\
\end{split}$ (104)

We have obtained all the surface and bulk terms needed to compute the final
expression for $n^{(6)}(\mathbf{k})$ (see Sects. 4.1 and 4.2). The explicit
calculation is very long and tedious (it needs substitution of eqs (82)). 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

$\begin{split}&\sum_{\gamma=\perp,\eta}\Big{(}\mathfrak{S}_{\gamma}^{(3)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)}(\mathbf{k})\Big{)}\Big{(}\mathfrak{S}_{\gamma}^{(3)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(3)\ast}(\mathbf{k})\Big{)}\\\
=&\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4}}\mathcal{H}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4})f^{ab_{1}b_{2}}f^{db_{3}b_{4}}\rho_{P}^{b_{1}}(\mathbf{p}-\mathbf{p}_{2})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{q}-\mathbf{p}_{4})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{ae}(\mathbf{k}-\mathbf{p})U^{de}(-\mathbf{k}-\mathbf{q})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{4}}\mathcal{H}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})f^{db_{3}b_{4}}f^{ebc}\rho^{b_{1}}_{P}(\mathbf{p}_{1})\rho^{b_{2}}_{P}(\mathbf{p}_{2})\rho^{b_{3}}_{P}(\mathbf{q}-\mathbf{p}_{4})\rho^{b_{4}}_{P}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{b_{1}b}(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})U^{b_{2}c}(\mathbf{p}-\mathbf{p}_{2})U^{de}(-\mathbf{k}-\mathbf{q})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}}\mathcal{H}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})f^{abc}f^{ade}\rho_{P}^{b_{1}}(\mathbf{p}_{1})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{p}_{3})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{b_{1}b}(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})U^{b_{2}c}(\mathbf{p}-\mathbf{p}_{2})U^{b_{3}d}(-\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})U^{b_{4}e}(\mathbf{q}-\mathbf{p}_{4})\\\
&+{\rm c.c.}\,.\end{split}$ (105)

The three scalar functions
$\mathcal{H}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4})$,
$\mathcal{H}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})$,
$\mathcal{H}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})$
only depend on transverse momentum. We present their explicit expressions in
the following.

$\mathcal{H}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4})=-\frac{1}{4}\frac{(\mathbf{p}_{2}\times\mathbf{p})}{|\mathbf{p}-\mathbf{p}_{2}|^{2}p_{2}^{2}}\frac{(\mathbf{p}\cdot\mathbf{q})}{p_{\perp}^{2}q_{\perp}^{2}}\frac{(\mathbf{p}_{4}\times\mathbf{q})}{|\mathbf{q}-\mathbf{p}_{4}|^{2}p_{4}^{2}}.$
(106)

We have used the identity
$(\mathbf{p}_{2}\times\mathbf{q})p_{\perp}^{2}-(\mathbf{p}\times\mathbf{q})(\mathbf{p}\cdot\mathbf{p}_{2})=(\mathbf{p}\cdot\mathbf{q})(\mathbf{p}_{2}\times\mathbf{p})$.
It is obvious that $\mathcal{H}_{1}$ is a real function of transverse
momentum. Up to a minus sign, $\mathcal{H}_{1}$ is nothing but the product of
the two effective vertex introduced in Eq. (98)

$\mathcal{H}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4})=-\Gamma^{j}_{1}(\mathbf{p}_{2},\mathbf{p})\Gamma^{j}_{1}(\mathbf{p}_{4},\mathbf{q})\,.$
(107)

See Fig. 6(a) for the momentum assignments.

There are six terms in $\mathcal{H}_{2}$,

$\begin{split}&\mathcal{H}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})\\\
=&\Bigg{\\{}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{p_{1}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{(\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}+\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\left(i\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}-1\right)\Bigg{\\}}\\\
&\quad\times\frac{1}{2}\frac{\mathbf{k}\times\mathbf{q}}{k_{\perp}^{2}q_{\perp}^{2}}\frac{\mathbf{p}_{4}\times\mathbf{q}}{|\mathbf{q}-\mathbf{p}_{4}|^{2}p_{4}^{2}}+\Bigg{\\{}-\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{p_{1}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\left[\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\cdot\mathbf{p}}{p_{\perp}^{2}}+\frac{1}{2}\frac{\mathbf{p}\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\right]\\\
&\quad+\frac{1}{2}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{(\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\left(i\frac{\mathbf{p}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}-1\right)\\\
&\quad-\frac{i}{2}\frac{(\mathbf{k}-\mathbf{p})\times\mathbf{p}_{1}}{p_{1}^{2}|\mathbf{k}-\mathbf{p}|^{2}}\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\frac{k_{\perp}^{2}+p_{\perp}^{2}-\mathbf{k}\cdot\mathbf{p}}{|\mathbf{k}\times\mathbf{p}|}\Bigg{\\}}\frac{1}{2}\frac{\mathbf{k}\cdot\mathbf{q}}{k_{\perp}^{2}q_{\perp}^{2}}\frac{\mathbf{p}_{4}\times\mathbf{q}}{|\mathbf{q}-\mathbf{p}_{4}|^{2}p_{4}^{2}}.\end{split}$
(108)

Using the effective vertices introduced in Eqs. (98) and (102),
$\mathcal{H}_{2}$ is precisely the product of these two vertices

$\mathcal{H}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})=\Gamma_{2}^{j}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p},\mathbf{k})\Gamma_{1}^{j}(\mathbf{p}_{4},\mathbf{q})\,.$
(109)

See Fig. 6(b) for the various momentum definitions. $\mathcal{H}_{2}$ contains
both real and imaginary parts. The imaginary part involves the single sign
function
$\mathrm{sgn}(\mathbf{k}\times\mathbf{q})=\mathbf{k}\times\mathbf{q}/|\mathbf{k}\times\mathbf{q}|$
and will be responsible for the asymmetry between $n^{(6)}(\mathbf{k})$ and
$n^{(6)}(-\mathbf{k})$.

There are thirteen terms in $\mathcal{H}_{3}$,

$\begin{split}&\mathcal{H}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})\\\
=&-\frac{1}{2k_{\perp}^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{(\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{k}+\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{q}-\mathbf{p}_{4})\cdot\mathbf{p}_{4}}{p_{4}^{2}}\\\
&+\frac{1}{2k_{\perp}^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{k}+\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}_{4}}{p_{4}^{2}}\\\
&-\frac{1}{2k_{\perp}^{2}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\left[\frac{\mathbf{k}\cdot\mathbf{p}}{p_{\perp}^{2}}-\frac{1}{4}\frac{(\mathbf{p}\times\mathbf{k})}{p_{\perp}^{2}}\right]\frac{(\mathbf{k}+\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\cdot\mathbf{p}_{4}}{p_{4}^{2}}\frac{(\mathbf{q}\times\mathbf{k})}{q_{\perp}^{2}}\\\
&-\frac{i}{k^{2}_{\perp}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{(\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{k}+\mathbf{q}+\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}_{4}}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(\frac{\mathbf{k}\cdot(\mathbf{k}+\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}-i\right)\\\
&-\frac{i}{2k^{2}_{\perp}}\frac{(\mathbf{k}-\mathbf{p})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\left[\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\cdot\mathbf{p}}{p_{\perp}^{2}}-\frac{1}{2}\frac{\mathbf{p}\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{(\mathbf{p}\times\mathbf{k})}{p_{\perp}^{2}}\right]\\\
&\qquad\times\Bigg{[}\frac{(\mathbf{k}+\mathbf{q}+\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{q}-\mathbf{p}_{4})\cdot\mathbf{p}_{4}}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(\frac{\mathbf{q}\cdot(\mathbf{q}+\mathbf{k})}{|\mathbf{q}\times\mathbf{k}|}+i\right)\\\
&\qquad\quad+\frac{(\mathbf{k}+\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}_{4}}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\frac{(k_{\perp}^{2}+q_{\perp}^{2}+\mathbf{k}\cdot\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}\Bigg{]}\\\
&-\frac{1}{2k^{2}_{\perp}}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{k}\times\mathbf{p}}{p_{\perp}^{2}}\frac{(\mathbf{k}+\mathbf{q}+\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}_{4}}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\\\
&\qquad\qquad\times\left(\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{p})}{|\mathbf{k}\times\mathbf{p}|}+i\right)\left(\frac{\mathbf{k}\cdot(\mathbf{k}+\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}-i\right)\\\
&-\frac{1}{4k^{2}_{\perp}}\frac{(\mathbf{k}-\mathbf{p})\times\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{\mathbf{p}\times\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{p}\times\mathbf{k}}{p_{\perp}^{2}}\frac{(\mathbf{k}+\mathbf{q}+\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{q}-\mathbf{p}_{4})\cdot\mathbf{p}_{4}}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\\\
&\qquad\times\frac{(k_{\perp}^{2}+p_{\perp}^{2}-\mathbf{k}\cdot\mathbf{p})}{|\mathbf{p}\times\mathbf{k}|}\left(\frac{\mathbf{q}\cdot(\mathbf{q}+\mathbf{k})}{|\mathbf{q}\times\mathbf{k}|}+i\right)\\\
&-\frac{1}{8k^{2}_{\perp}}\frac{(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{(\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{p_{2}^{2}}\frac{\mathbf{p}\times\mathbf{k}}{p_{\perp}^{2}}\frac{(\mathbf{k}+\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}_{4}}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\\\
&\qquad\qquad\times\frac{(k_{\perp}^{2}+q_{\perp}^{2}+\mathbf{k}\cdot\mathbf{q})}{|\mathbf{q}\times\mathbf{k}|}\left(\frac{\mathbf{p}\cdot(\mathbf{p}-\mathbf{k})}{|\mathbf{p}\times\mathbf{k}|}-i\right)\\\
&-\frac{1}{8k^{2}_{\perp}}\frac{(\mathbf{k}-\mathbf{p})\times\mathbf{p}_{1}}{|\mathbf{k}-\mathbf{p}|^{2}p_{1}^{2}}\frac{(\mathbf{p}\times\mathbf{p}_{2})}{p_{2}^{2}}\frac{\mathbf{p}\times\mathbf{k}}{p_{\perp}^{2}}\frac{(\mathbf{k}+\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}+\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{q}\times\mathbf{p}_{4})}{p_{4}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\\\
&\qquad\qquad\times\frac{(k_{\perp}^{2}+q_{\perp}^{2}+\mathbf{k}\cdot\mathbf{q})}{|\mathbf{q}\times\mathbf{k}|}\frac{(k_{\perp}^{2}+p_{\perp}^{2}-\mathbf{k}\cdot\mathbf{p})}{|\mathbf{p}\times\mathbf{k}|}.\\\
\end{split}$ (110)

(a)

(b)

(c)

Figure 6: Color structure for terms in the $|M_{(3)}|^{2}$ that contain two,
three and four Wilson lines. The dashed line separates the amplitude and the
complex conjugate amplitude. In (a) two gluon lines crossing the target
nucleus, (b) three gluon lines crossing the target nucleus and (c) four gluon
lines crossing the target nucleus.

It is an easy exercise to show that $\mathcal{H}_{3}$ can be expressed in
terms of the effective vertex we introduced in Eq. (102)

$\begin{split}&\mathcal{H}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})+\mathcal{H}_{3}^{\ast}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})\\\
=&-\Gamma_{2}^{j}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p},\mathbf{k})\Gamma_{2}^{j\ast}(\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{q},-\mathbf{k}).\\\
\end{split}$ (111)

See Fig. 6(c) for the various momentum assignments in the amplitude and the
complex conjugate amplitude.

There are two sign functions present in the expression of $\mathcal{H}_{3}$,
$\mathrm{sgn}(\mathbf{k}\times\mathbf{p})$ and
$\mathrm{sgn}(\mathbf{k}\times\mathbf{q})$. Imginary parts of
$\mathcal{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 $\mathbf{k}\leftrightarrow-\mathbf{k}$.

As for the color structure of the order-$g^{3}$ amplitude squared Eq. (105),
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. (98)
and (102), 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.

$\begin{split}&\sum_{\gamma=\perp,\eta}\mathfrak{S}_{\gamma}^{(1)\ast}(\mathbf{k})\Big{(}\mathfrak{S}_{\gamma}^{(5)}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(5)}(\mathbf{k})\Big{)}+\mathfrak{S}_{\gamma}^{(1)}(\mathbf{k})\Big{(}\mathfrak{S}_{\gamma}^{(5)\ast}(\mathbf{k})+\mathfrak{B}_{\gamma}^{(5)\ast}(\mathbf{k})\Big{)}\\\
=&\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4}}\mathcal{F}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})f^{dab_{3}}f^{ab_{1}b_{2}}\rho_{P}^{b_{1}}(\mathbf{p}-\mathbf{p}_{2})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{q}-\mathbf{p})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{de}(\mathbf{k}-\mathbf{q})U^{b_{4}e}(-\mathbf{k}-\mathbf{p}_{4})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}}\mathcal{F}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})f^{bb_{1}b_{2}}f^{ade}\rho_{P}^{b_{1}}(\mathbf{p}-\mathbf{p_{2}})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{p}_{3})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{ba}(\mathbf{q}-\mathbf{p})U^{b_{3}d}(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})U^{b_{4}e}(-\mathbf{k}-\mathbf{p}_{4})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}}\mathcal{F}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})f^{abc}f^{ade}\rho^{b_{1}}_{P}(\mathbf{p}_{1})\rho^{b_{2}}_{P}(\mathbf{p}_{2})\rho^{b_{3}}_{P}(\mathbf{p}_{3})\rho^{b_{4}}_{P}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{b_{1}b}(\mathbf{p}-\mathbf{p}_{1})U^{b_{2}c}(\mathbf{q}-\mathbf{p}-\mathbf{p}_{2})U^{b_{3}d}(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})U^{b_{4}e}(-\mathbf{k}-\mathbf{p}_{4})\\\
&+{\rm c.c.}\end{split}$ (112)

The three kinematic functions
$\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3}$ only depend on transverse
momentum. Their explicit expressions are presented in the following. For
$\mathcal{F}_{1}$,

$\mathcal{F}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})=\frac{1}{2q_{\perp}^{2}|\mathbf{p}-\mathbf{p}_{2}|^{2}|\mathbf{q}-\mathbf{p}|^{2}}\left(\frac{-(\mathbf{p}\times\mathbf{q})}{p_{\perp}^{2}}+\frac{\mathbf{p}_{2}\times\mathbf{q}}{3p_{2}^{2}}\right)\left(\frac{(\mathbf{k}\times\mathbf{q})}{k_{\perp}^{2}}+\frac{\mathbf{p}_{4}\times\mathbf{q}}{p_{4}^{2}}\right).$
(113)

Subtracting the Lipatov vertex in the complex conjugate amplitude

$L^{j}(\mathbf{p}_{4},-\mathbf{k})=\left(\frac{\mathbf{p}_{4}^{j}}{p_{4}^{2}}+\frac{\mathbf{k}^{j}}{k_{\perp}^{2}}\right),$
(114)

we introduce an effective vertex at order-$g^{5}$,

$\Upsilon^{j}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2})=\frac{1}{2|\mathbf{p}-\mathbf{p}_{2}|^{2}|\mathbf{q}-\mathbf{p}|^{2}}\left(\frac{\mathbf{p}\times\mathbf{q}}{p_{\perp}^{2}}-\frac{\mathbf{p}_{2}\times\mathbf{q}}{3p_{2}^{2}}\right)\frac{\epsilon^{ij}\mathbf{q}^{i}}{q_{\perp}^{2}}.$
(115)

Then the coefficient function can be expressed as a product of two effective
vertices

$\mathcal{F}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})=\Upsilon^{j}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2})L^{j}(\mathbf{p}_{4},-\mathbf{k})\,.$
(116)

See Fig. 7(a) for the momentum assignments. It is interesting to note that the
effective vertex $\Upsilon_{1}^{j}(\mathbf{p},\mathbf{q},\mathbf{q}_{2})$ is
independent of the momentum $\mathbf{k}$. It is also apparent that
$\mathcal{F}_{1}$ is a real function.

(a)

(b)

(c)

Figure 7: Color structures for the terms in $M_{(1)}^{\ast}M_{(5)}$ that
contain two, three and four Wilson lines. The dashed line separates the
amplitude and the complex conjugate amplitude. In (a) two gluon lines crossing
the target nucleus, (b) three gluon lines crossing the target nucleus and (c)
four gluon lines crossing the target nucleus.

There are eight terms in $\mathcal{F}_{2}$. The final result can also be
expressed as a product of two effective vertices

$\mathcal{F}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})=\Upsilon^{j}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})L^{j}(\mathbf{p}_{4},-\mathbf{k})$
(117)

We decompose the order-$g^{5}$ effective vertex $\Upsilon^{j}_{2}$ into two
parts

$\Upsilon^{j}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})=\Upsilon_{2}^{\|}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\frac{\mathbf{k}^{j}}{k_{\perp}^{2}}+\Upsilon_{2}^{\perp}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\frac{\epsilon^{ij}\mathbf{k}^{i}}{k_{\perp}^{2}}.$
(118)

The two functions are

$\begin{split}&\Upsilon_{2}^{\|}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\\\
=&-\frac{1}{2}\frac{\mathbf{p}_{2}\times\mathbf{p}}{|\mathbf{p}-\mathbf{p}_{2}|^{2}p_{2}^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}}\left(\frac{(\mathbf{k}-\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}-\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{q_{\perp}^{2}p_{3}^{2}}\right)\\\
&\qquad-\frac{1}{2}\frac{\mathbf{p}_{2}\times\mathbf{p}}{|\mathbf{p}-\mathbf{p}_{2}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(i\frac{\mathbf{k}\cdot\mathbf{q}}{|\mathbf{k}\times\mathbf{q}|}-1\right)\\\
&\qquad+\frac{1}{2}\frac{(\mathbf{p}_{2}\times\mathbf{p})}{p_{2}^{2}|\mathbf{p}-\mathbf{p}_{2}|^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{p}\cdot\mathbf{q})}{p_{\perp}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(i\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}-1\right)\\\
\end{split}$ (119)

and

$\begin{split}&\Upsilon_{2}^{\perp}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\\\
=&\frac{1}{2}\frac{(\mathbf{p}_{2}\times\mathbf{p})}{|\mathbf{p}-\mathbf{p}_{2}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{p}\cdot\mathbf{k})}{p_{\perp}^{2}}+\frac{1}{2}\frac{\mathbf{p}_{2}\times\mathbf{p}}{|\mathbf{p}-\mathbf{p}_{2}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})}{q_{\perp}^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}}\\\
&\qquad-\frac{1}{2}\frac{\mathbf{p}_{2}\times\mathbf{p}}{|\mathbf{p}-\mathbf{p}_{2}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(i\frac{\mathbf{q}\cdot(\mathbf{q}-\mathbf{k})}{|\mathbf{q}\times\mathbf{k}|}+1\right)\\\
&\qquad+i\frac{1}{2}\frac{(\mathbf{p}_{2}\times\mathbf{p})}{p_{2}^{2}|\mathbf{p}-\mathbf{p}_{2}|^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{p_{3}^{2}}\frac{(\mathbf{p}\cdot\mathbf{q})}{p_{\perp}^{2}}\frac{(\mathbf{k}\times\mathbf{q})}{q^{2}_{\perp}|\mathbf{k}-\mathbf{q}|^{2}}\frac{(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}.\end{split}$
(120)

See Fig. 7(a) for the momentum assignments. $\mathcal{F}_{2}$ contains both
real part and imaginary parts. The imaginary part is associated with the sign
function $\mathrm{sgn}(\mathbf{k}\times\mathbf{q})$.

There are nineteen terms in $\mathcal{F}_{3}$. Six of the terms involve the
auxiliary functions $\mathcal{I}_{1,2,3}$ and $\mathcal{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 $\mathcal{F}_{3}$
can also be expressed as a product of two effective vertices

$\mathcal{F}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})=\Upsilon^{j}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})L^{j}(\mathbf{p}_{4},-\mathbf{k})$
(121)

The third effective vertex at order-$g^{5}$, $\Upsilon^{j}_{3}$, can be
further decomposed into two parts depending whether it is parallel or
perpendicular to the momentum $\mathbf{k}$,

$\Upsilon^{j}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})=\Upsilon_{3}^{\|}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\frac{\mathbf{k}^{j}}{k_{\perp}^{2}}+\Upsilon_{3}^{\perp}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\frac{\epsilon^{ij}\mathbf{k}^{i}}{k_{\perp}^{2}}.$
(122)

Here the two functions $\Upsilon_{3}^{\|}$ and $\Upsilon_{3}^{\perp}$ are

$\begin{split}&\Upsilon_{3}^{\|}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\\\
=&\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{\perp}^{2}p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}i\mathcal{I}_{1}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
&+\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p^{2}_{\perp}p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}i\mathcal{I}_{2}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
&+\frac{\mathbf{p}\times\mathbf{p}_{1}}{p^{2}_{\perp}p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p})\times\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}i\mathcal{I}_{3}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
&+\left[-\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}}+\frac{1}{2}\frac{\mathbf{p}\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{q}\cdot\mathbf{p}}{p_{\perp}^{2}}\right]\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{q_{\perp}^{2}p_{3}^{2}}\\\
&-\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\left[\frac{1}{2}\frac{(\mathbf{k}-\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}+\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(i\frac{(\mathbf{k}\cdot\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}-1\right)\right]\\\
&+\left[-\frac{(\mathbf{p}\times\mathbf{p}_{1})}{p_{1}^{2}}\frac{\mathbf{q}\cdot\mathbf{p}}{p_{\perp}^{2}}+\frac{1}{2}\frac{\mathbf{p}\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{p}\times\mathbf{q})}{p_{\perp}^{2}}\right]\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(i\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}-1\right)\\\
\end{split}$ (123)

and

$\begin{split}&\Upsilon_{3}^{\perp}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{k})\\\
=&\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{\perp}^{2}p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{k}\times\mathbf{p}_{4})}{k^{2}_{\perp}p_{4}^{2}}i\mathcal{G}_{1}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
&+\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{\perp}^{2}p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p}-\mathbf{p}_{2})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{k}\times\mathbf{p}_{4})}{k^{2}_{\perp}p_{4}^{2}}i\mathcal{G}_{2}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
&+\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{\perp}^{2}p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p})\times\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{k}\times\mathbf{p}_{4})}{k^{2}_{\perp}p_{4}^{2}}i\mathcal{G}_{3}(\mathbf{p},\mathbf{q},\mathbf{k})\\\
&+\left[-\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{q}\times\mathbf{p}}{p_{\perp}^{2}}+\frac{1}{2}\frac{\mathbf{p}\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{q}\cdot\mathbf{p}}{p_{\perp}^{2}}\right]\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{k}\cdot(\mathbf{k}-\mathbf{q})}{q_{\perp}^{2}}\\\
&\quad+\left[-\frac{1}{2}\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{k}\cdot\mathbf{p}}{p_{\perp}^{2}}+\frac{1}{3}\frac{\mathbf{p}\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{p}\times\mathbf{k}}{p_{\perp}^{2}}\right]\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\\\
&\quad+i\left[-\frac{\mathbf{p}\times\mathbf{p}_{1}}{p_{1}^{2}}\frac{\mathbf{q}\cdot\mathbf{p}}{p_{\perp}^{2}}+\frac{1}{2}\frac{\mathbf{p}\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{p}\times\mathbf{q})}{p_{\perp}^{2}}\right]\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}\mathbf{p}_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q})\times\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{(\mathbf{k}\times\mathbf{q})}{q^{2}_{\perp}}\frac{(k_{\perp}^{2}+q_{\perp}^{2}-\mathbf{k}\cdot\mathbf{q})}{|\mathbf{k}\times\mathbf{q}|}\\\
&\quad-\frac{(\mathbf{p}-\mathbf{p}_{1})\cdot\mathbf{p}_{1}}{p_{1}^{2}}\frac{(\mathbf{q}-\mathbf{p})\cdot\mathbf{p}_{2}}{|\mathbf{q}-\mathbf{p}|^{2}p_{2}^{2}}\frac{(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})\cdot\mathbf{p}_{3}}{|\mathbf{k}-\mathbf{q}|^{2}p_{3}^{2}}\frac{\mathbf{k}\times\mathbf{q}}{q_{\perp}^{2}}\left(i\frac{\mathbf{q}\cdot(\mathbf{q}-\mathbf{k})}{|\mathbf{q}\times\mathbf{k}|}+1\right).\end{split}$
(124)

See Fig. 7(c) for the momentum assignments.

In addition to the explicit imaginary part that is associated with the sign
function $\mathrm{sgn}(\mathbf{k}\times\mathbf{q})$, the six functions
$\mathcal{I}_{1,2,3}$ and $\mathcal{G}_{1,2,3}$ are in general complex.
Therefore the effective vertex is complex.

Inspecting Eq. (112), we conclude that we only need to consider three terms as
far as the color structure is concerned. The color structures of the terms
containing two, three and four Wilson lines are shown in Fig. 7. Of course,
the complication is hidden in the three kinematic functions
$\mathcal{F}_{1},\mathcal{F}_{2},\mathcal{F}_{3}$. With the help of the
effective vertices, these coefficients functions can also be computed.

Adding Eq. (105) and Eq. (112), one obtains the main result of this paper, the
first saturation correction to single gluon production

$\begin{split}n^{(6)}(\mathbf{k})=&\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4}}\mathcal{H}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})f^{ab_{1}b_{2}}f^{cb_{3}b_{4}}\rho_{P}^{b_{1}}(\mathbf{p}-\mathbf{p}_{2})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{q}-\mathbf{p}_{4})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{ae}(\mathbf{k}-\mathbf{p})U^{ce}(-\mathbf{k}-\mathbf{q})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4}}\mathcal{F}_{1}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})f^{dab_{3}}f^{ab_{1}b_{2}}\rho_{P}^{b_{1}}(\mathbf{p}-\mathbf{p}_{2})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{q}-\mathbf{p})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{de}(\mathbf{k}-\mathbf{q})U^{b_{4}e}(-\mathbf{k}-\mathbf{p}_{4})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{4}}\mathcal{H}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{4},\mathbf{k})f^{db_{3}b_{4}}f^{ebc}\rho^{b_{1}}_{P}(\mathbf{p}_{1})\rho^{b_{2}}_{P}(\mathbf{p}_{2})\rho^{b_{3}}_{P}(\mathbf{q}-\mathbf{p}_{4})\rho^{b_{4}}_{P}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{b_{1}b}(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})U^{b_{2}c}(\mathbf{p}-\mathbf{p}_{2})U^{de}(-\mathbf{k}-\mathbf{q})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}}\mathcal{F}_{2}(\mathbf{p},\mathbf{q},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})f^{bb_{1}b_{2}}f^{ade}\rho_{P}^{b_{1}}(\mathbf{p}-\mathbf{p_{2}})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{p}_{3})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{ba}(\mathbf{q}-\mathbf{p})U^{b_{3}d}(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})U^{b_{4}e}(-\mathbf{k}-\mathbf{p}_{4})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}}\mathcal{H}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})f^{abc}f^{ade}\rho_{P}^{b_{1}}(\mathbf{p}_{1})\rho_{P}^{b_{2}}(\mathbf{p}_{2})\rho_{P}^{b_{3}}(\mathbf{p}_{3})\rho_{P}^{b_{4}}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{b_{1}b}(\mathbf{k}-\mathbf{p}-\mathbf{p}_{1})U^{b_{2}c}(\mathbf{p}-\mathbf{p}_{2})U^{b_{3}d}(-\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})U^{b_{4}e}(\mathbf{q}-\mathbf{p}_{4})\\\
&+\frac{1}{\pi}\int_{\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4}}\mathcal{F}_{3}(\mathbf{p},\mathbf{q},\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3},\mathbf{p}_{4},\mathbf{k})f^{abc}f^{ade}\rho^{b_{1}}_{P}(\mathbf{p}_{1})\rho^{b_{2}}_{P}(\mathbf{p}_{2})\rho^{b_{3}}_{P}(\mathbf{p}_{3})\rho^{b_{4}}_{P}(\mathbf{p}_{4})\\\
&\qquad\qquad\times
U^{b_{1}b}(\mathbf{p}-\mathbf{p}_{1})U^{b_{2}c}(\mathbf{q}-\mathbf{p}-\mathbf{p}_{2})U^{b_{3}d}(\mathbf{k}-\mathbf{q}-\mathbf{p}_{3})U^{b_{4}e}(-\mathbf{k}-\mathbf{p}_{4})\\\
&+c.c.\,.\end{split}$ (125)

There are a few comments about $n^{(6)}(\mathbf{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
$\rho_{R}^{b}(\mathbf{k})=\rho_{P}^{a}(\mathbf{p})U^{ab}(\mathbf{k}-\mathbf{p})$,
the two color structures are

$\begin{split}&-f^{abc}f^{ade}\rho_{R}^{b}(\mathbf{k}-\mathbf{p})\rho_{R}^{c}(\mathbf{p})\rho_{R}^{d}(-\mathbf{k}-\mathbf{q})\rho^{e}_{R}(\mathbf{q})=\Big{[}\rho_{R}(\mathbf{k}-\mathbf{p}),\rho_{R}(\mathbf{p})\Big{]}^{a}\Big{[}\rho_{R}(-\mathbf{k}-\mathbf{q}),\rho_{R}(\mathbf{q})\Big{]}^{a},\\\
&-f^{abc}f^{ade}\rho_{R}^{b}(\mathbf{p})\rho_{R}^{c}(\mathbf{q}-\mathbf{p})\rho_{R}^{d}(\mathbf{k}-\mathbf{q})\rho_{R}^{e}(-\mathbf{k})=\Big{[}\rho_{R}(\mathbf{p}),\rho_{R}(\mathbf{q}-\mathbf{p})\Big{]}^{a}\Big{[}\rho_{R}(\mathbf{k}-\mathbf{q}),\rho_{R}(-\mathbf{k})\Big{]}^{a}.\\\
\end{split}$ (126)

Here $\mathbf{k}$ is some given external momentum. The integration variables
are $\mathbf{p},\mathbf{q}$. It is obvious that one can not redefine
$\mathbf{p},\mathbf{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 $\int d^{2}\mathbf{k}$ in the game,
one can redefine $\mathbf{q}\rightarrow\mathbf{k}$ and
$\mathbf{k}\rightarrow-\mathbf{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)}(\mathbf{k})$ is expressed as a functional of the color charge
densities $\rho_{P}$ and $\rho_{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 $\rho_{P}$. However, the average over $\rho_{T}$ contains
four adjoint Wilson lines and is in general very complicated albeit possible
Blaizot:2004wv ; Dominguez:2012ad ; Shi:2017gcq . There also exists other
approximation schemes, see discussions in Refs. Lappi:2015vta ; Kovner:2017ssr
; Altinoluk:2018ogz . Event averaging will be studied in a follow-up paper.

As far as color structure is concerned, $n^{(6)}(\mathbf{k})$ can be simply
computed from the six diagrams in Fig. 6 and Fig. 7. These six diagrams
capture the correct dependence on $\rho_{P}$ and $\rho_{T}$. Of course, the
realistic coefficients $\mathcal{F}_{1,2,3}$ and $\mathcal{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)}(\mathbf{k})$
from the six diagrams .

Compared with the leading order gluon production $n^{(2)}(\mathbf{k})$, the
most important feature of the first saturation correction
$n^{(6)}(\mathbf{k})$ is that it includes nontrivial final state interactions.
This is explicitly shown in figs. 6(b), 6(c), 7(b) and 7(c). Gluons produced