First Saturation Correction in High Energy Proton-Nucleus Collisions: I. Time evolution of classical Yang-Mills fields beyond leading order

In high energy proton-nucleus collisions, the single- and double-inclusive soft gluon productions at the leading order have been calculated and phenomenologically studied in various approaches for many years. These studies do not take into account the saturation and multiple rescatterings in the field of the proton. The first saturation correction to these leading order results (the terms that are enhanced by the combination $\alpha_s^2 \mu^2$, where $\mu^2$ is the proton's color charge squared per unit transverse area) has not been completely derived despite recent attempts using a diagrammatic approach. This paper is the first in a series of papers towards analytically completing the first saturation correction to physical observables in high energy proton-nucleus collisions. Our approach is to analytically solve the classical Yang-Mills equations in the dilute-dense regime using the Color Glass Condensate effective theory and compute physical observables constructed from classical gluon fields. In the current paper, the Yang-Mills equations are solved perturbatively in the field of the dilute object (the proton). Next-to-leading order and next-to-next-to-leading order analytic solutions are explicitly constructed. A systematic way to obtain all higher order analytic solutions is outlined


Introduction
The two decades worth of experimental measurements at RHIC and, then, the LHC have provided many unexpected results, including strong evidence for the formation of a strongly coupled plasma of quarks and gluons in heavy-ion collisions at high energy. This plasma demonstrated properties of a nearly perfect fluid; this fact facilitated a theoretical description of the collision dynamics in the framework of hydrodynamics starting just about 1 fm/c after the heavy ion impact (see Refs. [1][2][3] and references therein).
The success of the hydrodynamic description, however, cannot be complete without a detailed understanding of the initial non-equilibrium state. The properties of this state go beyond the range of applicability of hydrodynamics but are crucial in fitting experimental data; the evolution of this state towards equilibrated thermal nearly perfect liquid have been one of the open theoretical problems being extensively studied [4][5][6][7][8]. One dominant mechanism describing the initial phase is based on the saturation framework [9][10][11][12], also widely known as the Color Glass Condensate (CGC). According to the framework, the high energy particle production and scattering processes are dominated by the classical gluon fields providing a background for systematic weak-coupling computation of quantum correction on top of it.
Under laboratory conditions, collisions of heavy-ions create probably the most optimal environment for probing quark-gluon plasma near equilibrium, but at the same time they are poorly suited to study the initial state particle production. This is because most of the observables in heavy-ion collisions are sensitive not only to initial state, but also to rather strong final sate interactions [13,14]. However, to uniquely map the transport properties of the plasma, it is critical to extract information about the initial state in collisions where the final state is better understood and the initial state is expected to play the dominant role. This necessitates probing a nucleus and a nucleon with the smallest projectiles: proton and ultimately electron. Theoretically, a controlled, first principle description of such asymmetric collisions (e-A, p-A, heavy-light nuclei) is not as complex as A-A collisions.
In the CGC framework, the key building block of soft gluon production in hadronic collisions is the single inclusive gluon cross section for a fixed configuration of the valence charges. Then the multi-gluon productions can be constructed from it iteratively. Analytical calculations of the single inclusive gluon production in asymmetric hadronic collisions at leading order in the color charge density of the dilute projectile (e.g. proton) have been done by various groups for more than two decades [15][16][17][18][19]. The leading order result takes into account the multiple rescatterings/saturation in the dense nucleus (target) to all orders while treats the collision partner (projectile) as dilute object. Beyond the leading order result, the first saturation corrections in the projectile to single-and double-inclusive gluons production were partially calculated in Refs. [20][21][22]. These incomplete results were sufficient to convincingly demonstrate that, in the CGC framework, the first saturation corrections are responsible for the generation of the odd azimuthal anisotropy [22,23] which was missing at the leading order.
In order to calculate the first saturation corrections to the single gluon production amplitude, the authors of Ref. [21] used the diagrammatic approach based on the light-cone perturbation theory and the eikonal approximation. Reference [21] provides a result for the order-g 3 single gluon production amplitude; the order-g 5 gluon production amplitude was not evaluated but it is needed to establish the complete first saturation corrections to the single inclusive gluon production in high energy proton-nucleus collisions. An alternative computational approach was adopted in Ref. [22]. The authors of Ref. [22] solve the classical Yang-Mills equations in the dilute-dense regime considering the projectile charge density as parametrically small. Particle production is then constructed from the classical gluon fields using the Lehmann-Symanzik-Zimmermann (LSZ) reduction formula. Although this approach was proven to be more powerful in helping organize calculations and extract the odd azimuthal anisotropy of double-inclusive gluon production, the authors of Ref. [22] also only computed order-g 3 production amplitude. The goal of this series of papers is to systematically complete the effort started in Ref. [22] and more specifically: 1) to solve for the next-to-leading order solutions of the classical Yang-Mills equations; 2) calculate the order-g 5 gluon production amplitude; 3) complete the first saturation corrections to the single-and double-inclusive gluon productions; 4) evaluate the early time-dependence of the energy-momentum tensor T µν (τ, x) and its correlation T µν (τ, x)T µ ν (τ , x ) .
Before proceeding with solving the classical Yang-Mills equations, we want to outline the role of the saturation corrections in high energy nuclear collisions. General perturbative corrections are terms expressed as power series expansions in the strong coupling constant g. In high energy nuclear collisions, the colliding hadrons are highly Lorentz contracted and the number of "valence" color sources per unit area as a random variable scales as N ⊥ ∼ √ A 1/3 with A the nuclear atomic number for a nucleus. For large A, at each order g n (n ≥ 1), there are terms that are enhanced by N ⊥ . The most enhanced term at each perturbative order is the saturation correction we aim to compute.
In order to illuminate the meaning of the saturation correction term, we will use smallx gluon distribution of a high energy proton as an example. The amplitude at order g and g 3 are schematically shown in Fig. 1. The order-g 5 amplitude is illustrated in Fig. 2. We start our discussions with the leading order. At order g, only one of the color source radiates a gluon, but there are N ⊥ color sources, thus the amplitude is proportional to gN ⊥ . The number of produced gluons is proportional to the amplitude squared which is parametrically g 2 N 2 ⊥ ∼ α s A 1/3 . At order g 3 , there are two types of diagrams, see Fig.  1. For the first kind, the gluon is radiated from one single color source; in this case, the amplitude is proportional to g 3 N ⊥ . For the second kind, the gluon originates from interaction of two color sources; in this case, the amplitude is proportional to g 3 N 2 ⊥ . The saturation correction only takes into account diagrams proportional to g 3 N 2 ⊥ , as they are most enhanced by the nuclear effects. For this term, the amplitude square is proportional to g 6 N 4 ⊥ ∼ α 3 s A 2/3 . Compared to the leading order contribution, it is higher order in α 2 s A 1/3 , which is what we have alluded to as the saturation correction. Now, at order g 5 , there are three types of diagrams, see Fig. 2. The first type of diagrams only involve one color source to radiate a gluon. The amplitude is proportional to g 5 N ⊥ . The second type comes with two color sources radiating gluons. Its amplitude is proportional to g 5 N 2 ⊥ . And, finally, the last type of diagrams with three color sources emitting gluons leads to the amplitude proportional to g 5 N 3 ⊥ . At order g 5 , the saturation correction only takes into Figure 1: Schematic diagrams showing perturbative corrections to small-x gluon distribution at order g and order g 3 . Saturation correction at order g 3 only takes into account the type of diagrams which are parametrically proportional to g 3 N 2 ⊥ .
⊥ + + Figure 2: Schematic diagrams showing perturbative corrections to small-x gluon distribution at order g 5 . Saturation correction at order g 5 only considers the types of diagrams proportional to g 5 N 3 ⊥ .
account the diagram that is proportional to g 5 N 3 ⊥ . For this term, the amplitude squared is proportional to g 10 N 6 ⊥ ∼ α 5 s A. Compared to the leading order term, it is parametrically higher order in (α 2 s A 1/3 ) 2 , which is the second order in terms of the saturation correction. The above discussions can only be formally applied to a large nucleus A 1 ensuring that the saturation corrections are leading compared to the other perturbative contributions. Superficially, in case of proton-nucleus collisions, the saturation corrections should not play any special role, since the nuclear atomic number for proton is A = 1. However, this is not completely right. At high energy, the number of color sources could still be large for at least two reasons. First, the proton wave-function can be in a rare configuration at the moment of collisions. The configurations like this are believed to be responsible for the high multiplicity events observed in high energy pA collisions in the experiments at the LHC and the RHIC. Second, the high energy evolution equations (BK and B-JIMWLK, see Refs. [24][25][26][27][28][29][30][31][32][33]) predict proliferation of the color charges; this ultimately leads to a universal high energy fixed point at which all hadrons look alike. To incorporate this general situation, it is more appropriate to reformulate the about counting in terms of the saturation scale Q s of the projectile instead of the nuclear atomic number.
In the CGC framework [34,35], the color sources responsible for gluon radiations are characterized by the random color charge density ρ a (x − , x), which represents "valence" partonic degrees of freedom. Specifically, in the McLerran-Venugopalan model, the color charge density is assumed to follow the Gaussian distribution with width µ(x − , x). The longitudinally integrated Gaussian width µ 2 (x) = dx − µ 2 (x − , x) has the physical meaning of color charge squared per unit transverse area. It represents the Gaussian width of the random variable ρ a (x) = dx − ρ a (x − , x). As a random variable, the characteristic scale of the color charge density is ρ a (x) ∼ √ A 1/3 . On the other hand, the saturation scale is shown to be related to µ 2 by Q 2 s ∼ (g 2 µ) 2 ∼ α 2 s A 1/3 [36]. This highlights the fact that the saturation corrections represent an expansion in terms of the projectile's saturation scale squared. It is not surprising that the power counting is captured exactly by solving the Yang-Mills equations for the classical gluon fields; this was explicitly demonstrated by Kovchegov [37].
To sum up, the single inclusive soft gluon productions in high energy nuclear collisions can be parameterized as a double Taylor series expansion of the saturation scales of the projectile and the target, Q s,P and Q s,T , respectively [21,23,38]. Using the notation of Ref. [21], In the case of a dilute projectile and a dense target, resummation over the target saturation corrections is possible and Eq. (1.1) can be written as an expansion in the projectile saturation momentum The leading order result f 1 (Q 2 s,T /k 2 ⊥ ) is known, see Refs. [15,18]. Corrections f 2 and f i>2 are not known analytically at present.
In case of double-inclusive gluon production, we have, schematically Here for simplicity we consider |k 1 | = |k 2 | = k ⊥ . The leading order result, h 1 , was derived in Refs. [39,40]. The first saturation correction was computed partially -only the odd component under the transformation k 1 → −k 1 was extracted in Ref. [22,23]. The goal of this series of papers is to compute the complete first saturation corrections, that is f 2 (Q 2 s,T /k 2 ⊥ ) and h 2 (Q 2 s,T /k 2 ⊥ ). There are several reasons why saturation corrections are important. On the practical side, leading order result does not include any final state interactions of the produced gluons. It basically assume the gluons propagate freely once created at proper time τ = 0. This might be a reasonable assumption for a dilute system created. The first saturation correction introduces non-trivial gluon interactions through three-gluon and four-gluon vertices. These interactions might be responsible for the onset of isotropization, thermalization and hydrodynamization. Additionally, having an expression for the first saturation corrections provides direct comparison of the relative importance of the initial state vs. final state effects. Furthermore, as alluded before, the most important feature of the final state effects is the generation of odd harmonics in multigluon distributions. Finally, the first saturation correction allows one to estimate the role of higher order contributions in the dilute-dense approximation. In all these cases, the saturation corrections are indispensable for any attempts to compare theory with experimental data.
On the academic side, rigorously calculating the first saturation correction is a first step towards including all saturation corrections. The ultimate goal is to resum all order saturation corrections and thus solve the dense-dense scattering problem analytically [41][42][43]. This is one of the unsolved problems in high energy QCD.
It should be mentioned that the classical Yang-Mills equations can be and were solved numerically. This approach was used for calculating the single-and double-inclusive gluon productions to all orders in both projectile and target color charge densities [44][45][46]. These calculations however rely on truncating the final state interactions at proper finite time.
The paper is organized as follows. After a brief introduction of the CGC framework and the classical Yang-Mills equations in sec. 2, we discuss the initial conditions in sec. 3. This includes explicit derivation of high order expansions; we also discuss a few convenient forms of different gauge fixings. The subsequent sections solve the classical Yang-Mills equations at orders g, g 3 and g 5 using the method of variation of parameters and Garf's formula for Bessel functions. For the discussions in sec. 8, we review which phsyical quantities can be obtained using our results.

The Color Glass Condensate Effective Theory
The Color Glass Condensate (CGC) effective theory concerns quantum chromodynamics in the hight energy limit [47][48][49]. It is based on a formal separation of large and small longitudinal momentum modes of partons inside a hadron. The partonic degrees of freedom with large longitudinal momenta (large-x) are effectively described by the color charge density ρ a (x). The gluons with small longitudinal momentum (small-x) are the dominant degrees of freedom and they are characterized by the classical gluon fields A µ (x). The color charge density is responsible for the production of the gluon fields through the classical Yang-Mills equations with the covariant derivative D µ = ∂ µ − igA µ and the field strength tensor For a right-moving hadron at high energy, the color current J µ (x) = δ µ+ ρ(x − , x) is approximately independent of light-cone time x + as far as the dynamics of the small-x gluons is concerned.
In applying CGC to high energy nuclear collisions with a right-moving projectile and a left-moving target, the color current can be approximated as J µ = δ µ+ ρ P (x − , x) + δ µ− ρ T (x + , x). Before the collisions, two sheets of small-x gluons, generated separately by the projectile and the target, approach each other at the speed of light, see Fig. 3. The collisions happen instantaneously (high-energy approximation). After the collisions, the large-x color charges are still approximately traveling along the lightcone x + = 0 and x − = 0 while classical gluon fields are produced in the forward lightcone x + > 0, x − > 0. The dynamics of the produced gluon fields is governed by the sourceless Yang-Mills equations.
A µ (x) =? The initial conditions are crucial as they encode the information about the instantaneous collisions. Once this initial value problem for the classical Yang-Mills equations is solved, one can compute physical observables that depend on classical gluon fields. Eventually, through the initial conditions, any physical observable will be a functional of the color charge densities of the projectile and the target O(ρ P , ρ T ).
The event/initial configuration-averaged results are obtained by evaluating the average over projectile and target color charge densities separately. In the McLerran-Venugopalan model, the color charge densities are assumed to follow Gaussian distributions. Their two-point correlation functions are To solve the sourceless classical Yang-Mills equation in the forward light cone x + > 0, x − > 0, we follow the literatures [18,50] and consider the Fock-Schwinger gauge x − A + + x + A − = 0. In this gauge, the solutions can be parameterized as Note that boost-invariance is assumed so that the classical gluon fields are independent of the rapidity η. The coordinate system used is denoted by (τ, η, x) with τ = √ 2x + x − and η = 1 2 ln x + x − . In terms of α(τ, x), α i (τ, x), the classical Yang-Mills equations become The initial conditions were derived in Refs. [50,51]: (2.5) Here α i P (x) and α i T (x) are the Weizsacker-Williams gluon fields of the projectile and target, respectively. The fields α i P (x) and α i T (x) are two dimensional pure gauge fields; they depend on the transverse coordinate and can be parameterized using Wilson lines with V (x) = exp{igΦ P (x)} and U (x) = exp{igΦ T (x)}. The relation between Φ P (T ) and ρ P (T ) will be explained in more details in the following sections. Eqs. (2.4), (2.5) define an initial value problem for a set of second order partial differential equations. The goal of the paper is to solve the Yang-Mills equations in the dilute-dense regime relevant to high energy proton-nucleus collisions. In the case of protonnucleus scatterings, the color charge density of the proton is parametrically small ρ P ∼ g while the color charge density of the nucleus is dense ρ T ∼ 1/g. We proceed to solve the classical Yang-Mills equations by expressing the gluon fields as power series expansions in terms of ρ P . The Yang-Mills equations can then be solved order by order perturbatively. The dependence on ρ T is resummed to all orders through the Wilson line. In the next section, to be consistent with the expansions of Yang-Mills equations, the initial conditions will also be expressed as power series expansions in ρ P . It is worth pointing out that the expansion here is not the same as the conventional expansion in the strong coupling constant g, as we only keep terms that are enhanced by ρ P at each order of perturbative expansion in g.

Expanding the Weizsacker-Williams field
From Eq. (2.6), the field generated by the "valence" color charge densities of the proton, also known as the Weizsacker-Williams (WW) field, can be expanded as (3.1) To simplify the notation, we dropped the subscript of Φ(x) for the projectile. The structure of the equation suggests that the WW gluon field in Eq. (3.1) is a pure gauge field. Additionally, it has to satisfy the static Yang-Mills equation Note that we have explicitly separated the g dependence in ρ P and thus parametrically ρ P ∼ 1 should be understood in the following discussions. This equation of motion constrains Φ(x), We solve for α i P (x) and Φ(x) perturbatively in g. Only terms with odd powers of coupling constant are nonvanishing Substituting the expansion of Φ(x) into Eq. (3.3), we obtain (3.5) To economize notation we introduced φ = 1 ∂ 2 ρ P . We only need up to order-g 5 expansions for the purpose of calculating the first saturation correction to gluon productions in high energy proton-nucleus collisions. Now it is straightforward to obtain the perturbative expressions for the WW gluon field order by order. Indeed, substituting Eq. (3.5) into Eq. (3.1), we get at the leading order. From this we conclude that the right hand-side of Eq. (3.2) is saturated automatically at the leading order, i.e. ∂ i α (1),i P = ρ P . Therefore all higher order contributions α (n),i P with (n ≥ 3) have to have vanishing gradient: We use this condition to cross-check the trivial algebra when deriving the higher order terms in the expansion of the WW gluon field. The cubic and quintic orders are (3.9) We factorized out the projection operator δ ij −∂ i ∂ j /∂ 2 to make the property (3.7) manifest.
We illustrate the expansion using Feynmann diagrams corresponding to terms at each order of the expansion. This is useful to ultimately establish a comparison with results from the diagramatic approach. Previously, this has been explicitly done in Ref. [37]; however, the discussion was limited to two color sources. Additionally, the expansion in Ref. [37] was conducted in the coupling constant g; in the current work, we perform the expansion in terms of the projectile saturation momentum, ∝ α s ρ P . Therefore our conclusions do not have to agree with Ref. [37] beyond order g 3 , as our goal is to account for the saturation correction rather than the perturbative correction! To proceed with the Feynmann diagrams it is convenient to explicitly define the inverse Laplacian operator appearing in φ( Here the scale Λ is an IR regularization scale. The vector potential is thus Now we are ready to proceed with the Feynmann diagrams. The leading order α i,(1) P is illustrated in Fig. 4. The order-g 3 WW gluon field α i,(3) P involves interactions of two color Figure 5: The order-g 3 WW gluon field α i,a,(3) P (x) produced at transverse position x by two color sources at transverse positions y 1 and y 2 .
charges. Written in explicit form, it has two terms (3.12) The Feynmann diagrams corresponding to these two terms are shown in Fig. 5. In these diagrams, the gluon emission at x from source at y 2 corresponds to the factor (x−y 2 ) i |x−y 2 | 2 while the factor ln (λ|x − y 1 |) illustrates one-gluon exchange between colored objects at x and y 1 .
The order-g 5 WW gluon field involves interactions of three color charges. In the expression for α i,(5) P , there are six topologically different contributions (3.13) We write them out explicitly one by one. The first term is (3.14) Figure 6: Two sets of diagrams contributing to order-g 5 WW field α i,a,(5) P (x) produced at transverse position x by three color sources at transverse position y 1 , y 2 and y 3 . Similar diagrams with different orderings of gluon exchanges are not shown.
The subsequent two terms are (3.16) These three terms correspond to the diagrams shown in Fig. 6a. Note that in the second diagram, there is an integration over all the possible transverse positions y.
The remaining three terms are The corresponding Feynman diagrams are shown in Fig. 6b.

Residual gauge fixing and initial conditions
The (3.20) There are a few choices of the sub gauges. We discuss them below.

Sub gauge transformation by U (x)
In the literature, when calculating particle production in pA collisions, U (x) was often chosen [18,22] to define the residual gauge fixing. With Ω(x) = U (x) the initial condition for transverse fields in Eq. (2.5) becomes Note that the target field α i T = i g U ∂ i U † is gauged away. In this form, both the gradient ∂ i ζ i and the curl ji ∂ j ζ i are nonvanishing. The initial condition for the longitudinal field in Eq. (2.5) becomes In this sub gauge, the order-g, order-g 3 and order-g 5 initial conditions are This sub gauge has the advantage that the initial conditions have clear physical meaning. The ζ i is the projectile WW gluon field α i P eikonally rotated by the target Wilson line U (x). The ζ is the difference between the eikonally rotated projectile WW gluon field and the gluon field generated by the eikonally rotated projectile color charge density. However, in solving the classical Yang-Mills equations for pA collisions beyond the leading order, this sub gauge is not the most convenient one. In a desirable sub gauge, either the gradient or the curl of ζ i (τ = 0, x) would vanish.

Sub gauge condition
It has been shown in Ref. [52] that the following sub gauge transformation guarantees that the lowest order gradient of β i vanishes, i.e. ∂ i β (1) i (τ = 0, x) = 0. To go beyond the lowest order, we consider the ansatz . Unitarity condition of Ω requires Σ † = Σ. Under this gauge transformation, the initial condition for the transverse fields in Eq. (2.5) becomes In obtaining the last equality, we used the Baker-Campbell-Hausdorff formula. We express Σ(x) as a power series expansion in terms of coupling constant g and solve Σ(x) order by order by imposing the requirement ∂ i β i (τ = 0, x) = 0. The results are (3.28) The initial conditions for the transverse field at the corresponding orders are (3.29) We will only need initial conditions up to order-g 5 . However, one can recursively obtain all higher order gauge transformations by imposing the condition ∂ i β i (n) (τ = 0, x) = 0 for n ≥ 7. It is not clear to us whether a closed form expression for Σ(x) exist or not.
On the other hand, the initial condition for the longitudinal field under the gauge transformation becomes From it, the order-g, order-g 3 , order-g 5 initial conditions for the longitudinal field are obtained (3.31) In obtaining the above results, we have used the fact that ∂ i α i P,(n) = 0 for n ≥ 3. The sub gauge condition ∂ i β i (τ = 0, x) = 0 resembles the general Coulomb gauge ∂ i β i (τ, x) = 0. Previously, in numerically solving the classical Yang-Mills equations, the Coulomb gauge condition was also used. However, instead of imposing it at τ = 0, it was imposed at some particularly chosen proper time τ 0 , at which physical observables were calculated [46,53].
One may wonder, whether it is possible to find a sub gauge transformation such that instead of the gradient of β i , the curl of β i is zero ij ∂ i β j (τ = 0, x) = 0. First of all, we want to point out that there is no Ω(x) that can completely gauge away The right hand side is a pure gauge field. On the other hand, both α i P and α i T are pure gauge fields, their sum cannot be a pure gauge field. Not even at the lowest order. It is not clear that the following condition on the curl of β i has a perturbative solution for Ω(x). As it will become clear in the following sections, the gradient of β i (τ, x) is an auxiliary field while its curl is a dynamical field. Therefore, as far as gluon production is concerned, the initial time Coulomb gauge constraining the gradient of β i (τ = 0, x) serves as a convenient gauge choice.
Our motivation for choosing different sub gauge transformations was purely to simplify computations as physical observables are gauge invariant. However, when discussing the time evolution of the gluon field (a gauge variant object), the concept of initial vs final state effects becomes blurred and, strictly speaking, not well-defined. A sub-gauge transformation may (and does) shift some final state effects to the realm of initial state effects and vice-versa. As a matter of fact our motivation was exactly to simplify the time evolution and transfer the computational burden to the initial conditions.
In the main body of this paper, the classical Yang-Mills equations are solved in the initial time Coulomb sub gauge ∂ i β i (τ = 0, x) = 0. In the Appendix D, solutions in the sub gauge determined by U (x) are given. The two sub gauges are related by W(x). By comparing gluon fields in the two gauges, it will become clear how final state interactions in one gauge become initial state effects in another gauge.

The Dynamical Equations and the Constraint Equation
As discussed in the previous section, the classical Yang-Mills equations are invariant under the sub gauge transformation. Thus the equations of motion in the initial time Coulomb sub gauge are obtainted by simple replacements of α, α i with β, β i in Eqs. (2.4).
We have written out the detailed expressions for the equations using the covariant derivative We also separated the linear and nonlinear terms in the equations and introduceβ = τ β to simplify the notation. The second equation in Eqs. (4.1) is first order in time derivative and thus serves as a constraint equation. Only two of the three field components β(τ, x), β i=1,2 (τ, x) are independent. To explicitly demonstrate this, it is convenient to split the transverse field into the gradient part and the curl part Here il is the two dimensional Levi-Civita symbol. Using ih ∂ h β i = ∂ 2 χ and ∂ i β i = ∂ 2 Λ, one can separate the third equation in Eqs. (4.1) into two equations: The constraint equation only imposes restriction on Λ The independent degrees of freedom areβ(τ, x) and χ(τ, x). The Λ(τ, x) is a non-dynamical field. In the Appendix A, it is proved perturbatively that the second order equation (4.3) is just a consequence of the first order equation (4.5). Thus, although the superficial appearance does not suggest it, the second order differential equation for Λ is not an independent dynamical equation.
In the following sections, we seek solutions of the classical Yang-Mills equations in terms of power series expansion in the coupling constant g.

Order-g Solutions
The order-g classical Yang-Mills equations are Performing the decomposition β , we trivially establish that the constraint equation together with the initial condition ∂ i β (1) i (τ = 0, x) = 0 is equivalent to Λ (1) (τ, x) = 0, as expected at this order.
The Yang-Mills equations are easier to solve in transverse momentum space. Using the convention for the Fourier transformations, we obtain that after the projection ih ∂ h β (1) i = ∂ 2 χ (1) the two independent equations are Here s ≡ k ⊥ τ and k ⊥ = |k| is the magnitude of the two dimensional transverse momentum. These two equations are easily recognized as standard Bessel equations. We require the solutions to be finite at τ = 0, only Bessel functions of first kind satisfy this constraint. We thus have where b η,⊥ (k) are fixed by the initial conditions j (τ = 0, k).

(5.6)
It is easy to check that ∂ i β (1) i (τ, x) = 0 for any τ . The order-g solution has already been obtained in Ref. [18]. The order-g equations in this sub-gauge are free field equations and the solutions are free field solutions. Somewhat nontrivial time dependence characterized by J 1 (k ⊥ τ )/k ⊥ τ and J 0 (k ⊥ τ ) is solely due to the Milne coordinates (τ, η, x). Now, we turn to higher orders and for the first time we will obtain order-g 3 and order-g 5 solutions.
6 Order-g 3 Solutions The order-g 3 classical Yang-Mills equations are with the initial conditions The first equation when transformed into momentums space is an inhomogeneous Bessel equation (s = k ⊥ τ ) In order to solve the inhomogeneous differential equations as Eq. (6.4), one can apply a well established method which is often referred to as variation of parameters. The method is briefly reviewed in the Appendix B. The two independent solutions for the corresponding homogeneous Bessel equation are J 1 (x) and Y 1 (x), whose Wronskian is W (x) = 2 πx . So the general solutions for the inhomogeneous equation can be formally expressed as The initial condition is finite at τ = 0, thus Y 1 (x) does not contribute, i.e. C 2 = 0. The coefficient C 1 can be then determined straightforwardly Putting everything together we get Further evaluations of the formal solution Eq. (6.6) require computing time integrals involving products of three Bessel functions in the integrand Here we face a difficulty because for integrands involving products of three or more Bessel functions, there are no known formula to compute the indefinite integrals. This is in a stark contrast to the case with integrands of only two Bessel functions:  This defines our strategy: we will aim at reducing the number of Bessel functions in the integrand from three (or more) to two in order to use the above equations to evaluate time integrals. The key step is to expressing a product of two Bessel functions in terms of an integral of one Bessel function using Graf's formula. Mathematical details are given in the Appendix C.
Substituting the expression of S η into the formal solution, one obtains (6.11) Introducing the notation w 1 = p ⊥ /k ⊥ and w 2 = |k − p|/k ⊥ , and using the formula and e iΨ = w 2 − w 1 cos φ w + i w 1 sin φ w (6.14) we can express the integrals of products of three Bessel functions as integrals of products of two Bessel functions. We have to pay a price of introducing an auxiliary angular integral. Nevertheless, after performing this manipulation, the structure of the solution simplifies significantly. Finally we will end up with an expression of the following form where we performed further simplification by using that fact that the Wronskian is πs . Note that w = 1 is a removable singularity; indeed, 2 π which is finite and well-defined.
Collecting everything together we obtain the final solution where The are a few notable features of the solution.
• First is that the time-dependent factors are completely determined by one type of Bessel function; in this case, it is Bessel function of first kind of order one J 1 (λτ )/λτ . However, the Bessel function contributes with different arguments. For the first term, the argument of the Bessel function is completely determined by the external momentum k. The second term is more involved, for given momenta k and p, the time-dependent factor sums over all possible momentum mode between w max = p ⊥ + |k − p| and w min = |p ⊥ − |k − p||. It should be pointed out that the second term has a removable singularity at w ⊥ = k ⊥ . It can be checked by performing Taylor expansions of the difference J 1 (w ⊥ τ )/w ⊥ τ − J 1 (k ⊥ τ )/k ⊥ τ ; it starts from w 2 ⊥ − k 2 ⊥ . This combination J 1 (w ⊥ τ )/w ⊥ τ − J 1 (k ⊥ τ )/k ⊥ τ also guarantees that the second term is zero at τ = 0.
Naively, it is expected that at asymptotically large τ , the gluon system should behave like free gas of gluons. This is not that easy to confirm on the level of the field; the asymptotic behavior of the Bessel function reads and thus the summation over all the momentum modes persists even at τ → ∞.
• The second feature is about the color structure of the solution. It involves interactions of two color charges in the proton. Let us look at each term in the solution in detail. First we can recognize from eq. (3.28) that is the order-g WW gluon field α i P,(1) eikonally color rotated by the target Wilson line U ad and then projected along the momentum k. Next, consider In the parenthesis, we have a difference between two terms. One is the order-g projectile WW gluon field eikonally rotated by the target Wilson line U ad . The other is the WW field generated by the eikonally rotated color density ρ a P U ad (it corresponds to the gluon cloud of the receding color charge or to the Fadeev-Kulish state). The net field is projected also along k for b η (k) and projected perpendicular to k for b ⊥ (k). In this sense, one can attribute b η (k) and b ⊥ (k) to two polarizations of the order-g gluon field produced in the collisions.
In the solution Eq. (6.17), the initial field at order-g 3 is Figure 7: Schematic representation of the color structure for the second term in β (3) (τ, k). The shaded bar represents the target nucleus. The eikonally rotated color charge density ρ R is represented using red lines.
The first term is the order-g 3 WW gluon field α i P, (3) which was eikonally color rotated by the target Wilson line U ac and then projected to the momentum k. The diagrams representing α i P, (3) have been shown in Fig. 5. The color structure of the second term in Eq. (6.17) can be schematically illustrated in Fig. 7.
We note that this discussion was specific for the used sub-gauge, since the color structure depends on sub gauge transformations. i (τ, k) As a first step, we want to demonstrate explicitly that from the order-g 3 constraint equation and the order-g Yang-Mills equations, the second order differential equation for ∂ i β is not a dynamical field. From the third equation of the set (6.1), the second order differential equation for ∂ i β

Solving for β
From the constraint equation of (6.1), one obtains In what follows, we will use the decomposition to separately solve for χ (3) (τ, x) and Λ (3) (τ, x).

(6.26)
The source term is given by In obtaining the last equality, we changed the integration variable p to k − p and used the fact that are antisymmetric under the exchange p ↔ k − p. Additionally, we applied the Bessel function identity J 2 (z) = 2J 1 (z)/z − J 0 (z), and we expressed product of two Bessel functions in terms of angular integral of one Bessel function: To obtain the solutions, directly integrating S Λ involves indefinite integrals of products of two Bessel functions of different orders and with different arguments. We are not aware of if these integrals can be done analytically. That is why we express products of two Bessel functions as an integral of one Bessel function.
The solution for Λ (3) (τ, k) is obtained by direct integration of S The time dependent factors are completely determined by Bessel function of first kind with order zero in the form 1 − J 0 (w ⊥ τ ). Again, for given k and p, all the momentum modes from |p ⊥ − |k − p|| to p ⊥ + |k − p| contribute to the argument of the Bessel function. Interestingly, the two polarization modes b ⊥ and b η do not mix.

The solution β
(3) The equation of motion for β (3) i (τ, k) is with the time dependent source term given by Instead of working with χ (3) (τ, k), it is convenient to project out the curl of β We changed variables p ↔ k − p, where appropriate, to symmetrize this expression. After the projection of Eq. (6.31), the equation of motion for the curl β This is an inhomogeneous differential equation and the corresponding homogeneous part is the the Bessel equation of the first kind. We again use the method of variation of parameters to solve it. The two independent general solutions for the corresponding homogeneous equation is J 0 (s) and Y 0 (s). Their Wronskian is W (s) = 2 πs . The formal solution of Eq. (6.35) is The solution is nonsingular at τ = 0, thus D 2 = 0. Substituting the explicit expression for S (3) ⊥ into the formal solution, one obtains We use the same strategy as described in the previous section for β (3) , i.e., we proceed by reducing integrals of three Bessel functions into integrals of two Bessel functions using Graf's formula. In this case, we have (6.38) Here w 1 = p ⊥ /k ⊥ , w 2 = |k − p|/k ⊥ and w 2 = w 2 1 + w 2 2 − 2w 1 w 2 cos φ. After performing this manipulations, the remaining integrals with two Bessel functions can be combined into the form Using these steps, the parts involving Bessel functions in Eq. (6.37) are simplified as 2 π 1 1 − w 2 (J 0 (ws) − J 0 (s)). (6.41) With this, we arrive at the final solution for β The time dependence is completely determined by one type of Bessel function -Bessel function of the first kind of zero order, J 0 (λτ ). The arguments of the Bessel functions are different, ranging from ||k − p| − p ⊥ | to |k − p| + p ⊥ .
For completeness, we supplement this solution with the detailed expression of the order-g 3 initial field (6.43)

Order-g 5 Solutions
The solutions presented in previous sections are sufficient to compute the full first saturation correction to single inclusive gluon production, as will be demonstrated in the second paper of this series [54]. However, in order to compute other interesting physical quantites like the energy-momentum tensor of the classical gluon fields, including the first saturation correction requires going beyond order-g 3 and finding order-g 5 solutions. Our motivation is to derive the energy-momentum tensor in a semi-analytic form to extract information about the energy density, pressure, stresses and initial flows after the collisions; these quantities can be as model initial conditions for a subsequent hydrodynamic evolution.
The order-g 5 equations of motion are and i ]]. Again, as in the previous section for g 3 -order, one can explicitly show that from the constraint equation for ∂ i β (5) i and all the lower order solutions, the second order differential equation for ∂ i β We want to write down the equations for the independent fields in momentum space. Performing the decomposition we can separate the equations for the curl and the gradient of β (5) i . In momentum space, we have The two dynamical equations governing the time evolutions of β (5) ⊥ (τ, k) and β (5) (τ, k) are i (τ, k).
The source terms in momentum space are and S We now have everything ready to find all the components of the gluon fields at order-g 5 .

Solving for β (5) (τ, k)
From eq. (7.9), the explicit expression of the source term S The last term involves three Bessel functions. The product of the three Bessel functions can be rewritten as To solve for β (5) (τ, k), we repeat the procedure used when solving for β (3) (τ, k). The formal solution obtained by the method of variation of parameters is with the coefficients C 1 and C 2 fixed by initial conditions at order-g 5 The time-dependent factors in each term of S (5) η (τ, k) can always be reduced to just one single Bessel function, Bessel function of first kind with order one J 1 (cτ ), although different terms might have different values of argument c. This reduction is done by using Graf's formula with c 2 = a 2 + b 2 − 2ab cos ψ. Next, the formula helps to carry out the integration over the source terms in the formal solution. Note that λ = c/k ⊥ . Schematically, our recipe is to replace each J 0 (aτ )J 1 (bτ ) factor in S To be specific, we need the following replacements (7.20) Using the above recipe, we obtained the order-g 5 solution β (5) (τ, k) The order-g 5 solution β (5) (τ, k) involves interactions of three color charges in the proton as evidenced by the double color commutators. From the definitions of b ⊥ (k), b η (k) in eqs. (6.21) and (6.20), the color structure of the three double commutators can be diagram- ] are shown in Figs. 8. As for the single commutators contain- R (x 1 ) + Figure 8: Schematic representation of the color structure for the double commutators in β (5) (τ, k). The shaded bar represents the target nucleus. The eikonally rotated color charge density ρ R is represented using red lines.
The time-dependent factors are expressed solely in terms of one type of Bessel functions J 1 (aτ )/aτ with the cost of introducing two auxilliary angular integrals. Typical terms also contain two transverse momentum integrations. These transverse momentum integrations reflect the momentum exchanges between the projectile and the target during the collisions. 7.2 Solving for β (5) ⊥ (τ, k) From eq. (7.10), the explicit expression of the source term S The product of three Bessel functions in the last two terms can be further reduced to a product of two Bessel functions using Graf's formula: To solve for β (5) ⊥ (τ, k), we again follow the same procedure used when solving for β (3) ⊥ (τ, k). The formal solution is obtained by using the method of variation of parameters β (5) with the coefficients determined by the initial condtions at order-g 5 The time-dependent factors in S ⊥ (τ, k) can always be reduced to one single type of Bessel function, Bessel function of the first kind with order zero J 0 (cτ ), although different terms might have different values of the argument c. To be more precise, using Graf's formula, only two possibilities are involved J 0 (cτ ), cos ψJ 0 (cτ ) (7.26) with c 2 = a 2 + b 2 − 2ab cos ψ. The next step is to use the formula (λ = c/k ⊥ ) It is clearly by now that the recipe is to replace time-dependent factors J 0 (aτ )J 0 (bτ ) and J 1 (aτ )J 1 (bτ ) in the source term S To be specific, we only need the following replacements in the source term The definitions of w ⊥ , u ⊥ , φ , θ are given in Eqs. (7.20).
The final expression for the order-g 5 solution β It is apparent that the transverse gluon field at order-g 5 solely depends on one type of Bessel function J 0 (aτ ) although different terms have different arguments. Color structure of the solution can be similarly analyzed as having done for β (5) (τ, k).

Solving for Λ (5) (τ, k)
Unlike the solutions β (5) (τ, k) and β (5) ⊥ (τ, k) , which are determined through the method of variation of parameters, the non-dynamical field Λ (5) (τ, k) is obtained by direct integration over the source term S (5) Λ (τ, k). In the following, we reorganize the time-dependent factors in S (5) Λ (τ, k) so that they only involve one type of Bessel function, Bessel function of the first kind with order one J 1 (cτ ). The expression of S (5) Λ (τ, k) in Eq. (7.11) can be further expressed as We compute each term separately. The first and the third terms can be combined together (7.32) Here We have also used the identity π −π dφ 2π The second and the fourth terms can be combined together Finally the fifth and the sixth terms are combined We used relation J 2 (z) = 2 z J 1 (z) − J 0 (z) to obtain (7.37) Using the explicit expression of S Λ (τ, k), the solution Λ (5) (τ, k) is ) .

(7.38)
Let us summarize and comment on the general procedures for solving the classical Yang-Mills equations perturbatively in the dilute-dense regime. At each fixed order g 2m+1 , the dynamical fields β (2m+1) (τ, k) and β (2m+1) ⊥ (τ, k) satisfy the inhomogeneous Bessel differential equations of orders one and zero, respectively. They are solved using the wellestablished method of variation of parameters. The success of this method relies on recombining the time-dependent factors in the source terms S only depends on J 0 (cτ ). Owing to the Graf' s formula, this is always possible. As for Λ (2m+1) (τ, k), it satisfies first order differential equation with source term S (2m+1) Λ . All one needs to do is express the time-dependent factors in S (2m+1) Λ (τ, k) in terms of Bessel functions J 1 (cτ ). Each time Graf's formula is used, an extra angular integral is introduced. Unfortunately, the number of terms in the solutions increases dramatically as one goes to higher and higher perturbative orders and the problem quickly becomes unmanageable analytically.

Discussions and Outlooks
In this paper, we have presented the first step towards completing the calculations of the first saturation corrections to physical observables in high energy proton-nucleus collisions.
We solved the classical Yang-Mills equations in the dilute-dense regime beyond the leading order. We explicitly constructed the order-g 3 and order-g 5 solutions. The main results are presented in Eqs. (6.17), (6.30), (6.42), (7.21), (7.30) and (7.38). The major mathematical technique that makes the analytic solutions possible is Graf's formula, which expresses product of two Bessel functions in terms of an angular integral of one Bessel function. As a consistence check of our main results, when the target Wilson line U (x) = 1, the gluon fields vanish β(τ, k) = 0 and β i (τ, k) = 0, i.e. there is no gluon production in the absence of scattering as expected.
There are a few apparent features of the solutions. First of all, the time dependent factors in the longitudinal gluon field are completely determined by J 0 (aτ ) again with possible different arguments a. Second, the order-g solutions β (1) (τ, k), β i (τ, k) do not involve mutual interactions of the "valence" color charges in the proton. Each color charge in the proton independently scatter on the target. On the other hand, the order-g 3 solutions β (3) (τ, k), β i (τ, k) represent interactions of two color charges in the proton. Their interactions can happen before or after the collisions with the target. The order-g 5 solutions β (5) (τ, k), β (5) i (τ, k) represent interactions of three color charges in the proton. Their interactions likewise can happen before or after the collisions with the target. It can also be that two color charges interact with each other before scattering on the target and then interact with the third color charge only after the collisions with the target.
Another important feature is related to the gauge dependence of the concepts of initial state effects and final state effects as defined on the basis of the field. In the main context of the paper, the solutions are given in the initial time Coulomb sub gauge. In appendix D, we present the order-g and order-g 3 solutions in the non-Coulomb sub gauge. These two sets of solutions are related by a gauge transformation. Performing direct comparison of these two sets of solutions, would convince the reader that some final state effects in the non-Coulomb subgauge become the initial state effects in the Coulomb subgauge.
On the other hand, physical observables are independent of gauge transformations. With the solutions of the classical Yang-Mills equations at hand, one can calculate several interesting physical quantities that can be constructed from the classical gluon fields. For example, the the energy-momentum tensor of the gluon fields produced in high energy pA collisions by T µν (τ, x) = F µλ (τ, x)F ν λ (τ, x) + 1 4 g µν F κλ (τ, x)F κλ (τ, x). Tracing over the color matrix is understood in the above definition like AB = A a B a = Tr(AB) with a = 1, . . . N 2 c − 1.
In the second paper of this series, we will calculate the single inclusive soft gluon production and double inclusive soft gluon production. These observables are constructed using the LSZ formulaâ a † p (τ ) = −iτ and taking the limit τ → ∞. Here H 0 (x) and H 1 (x) are Hankel functions of second kind, order zero and order one, respectively. The left right derivative is defined as . Unlike the energy-momentum tensor which is gague invariant by construction, the single inclusive gluon production byâ a, † pâ a p +ĉ a, † pĉ a p are explicitly dependent on β(τ, k) and β ⊥ (τ, k) which are gauge variant objects. Thus gauge invariance of the single inclusive production will serve as a non-trivial of the derived solutions.
It should be mentioned that to solve the classical Yang-Mills equations in the densedense regime, another semi-analytical approach was developed in Refs. [55,56]. In this method, the solutions are expanded as power series expansions in proper time τ . Although at each order in τ , all the saturation effects are included, the solutions are only meaningful when the values of τ are small. This small-τ expansion method cannot be used to rigorously calculate the single inclusive gluon production by the LSZ formula which requires taking the τ → ∞ limit. Our expansions in terms of coupling constant g are valid for all the proper time but can only take into account the saturation effects order by order. Amusingly one can combine these two analytical approaches into a double expansion in τ and g which gives a non-trivial insights into the dense-dense regime. We defer further discussion for a future publication. derived. Taking time derivative of Eq. (A.2) w.r.t. τ , one obtains For a general second order inhomogeneous differential equation a(x)y + b(x)y + c(x)y = f (x), (B.1) there are four steps to obtain its general solutions: • find two independent solutions y 1 (x) and y 2 (x) of the homogeneous equation The general solution for the homogeneous equation is where C 1 snd C 2 are coefficients to be fixed by the initial (and/or boundary ) conditions.
f (z) a(z) dz. (B.5) • finally, by adding together the general solution of the homogeneous equation y h (x) and the particular solution of the inhomogeneous equation y p (x), obtain the general solution for the inhomogeneous equation y(x) = C 1 y 1 (x) + C 2 y 2 (x) + x y 1 (z)y 2 (x) − y 1 (x)y 2 (z) C Integrals for the products of two Bessel functions There are a few general identities expressing the product of two Bessel functions in terms of an integral of one Bessel function from Dixion and Farrar's paper in 1933 [57]. For our problem here, only integer orders of Bessel functions are involved. The more well-known Graf's formula [58,59] serves our purpose: In this case, w = b 2 + a 2 − 2ab cos φ is symmetric with respect to a, b. and Here are a few examples that were used in the main body of the paper: It can be easily checked that φ + π + Ψ = Ψ . Thus the two integral representations are equivalent. Additionally, since w is an even function of φ, only φ−even part of the exponentials contributes: J 0 (az)J 1 (bz) =2 In the second equality, we have changed the integration variable from φ to w using cos φ = a 2 + b 2 − w 2 2ab , (C.13) dφ = wdw ab sin φ = 2wdw (C.14) The case of n = −2, m = 1 will yield the same expression. We can further simplify to get dφ 2π e −iφ e iΨ J 1 (wx), (C.23) but they are equivalent due to cos (φ − Ψ ) = − cos (2φ + Ψ).