Time-dependent observables in heavy ion collisions. Part I. Setting up the formalism

We adapt the Schwinger-Keldysh formalism to study heavy-ion collisions in perturbative QCD. Employing the formalism, we calculate the two-point gluon correlation function G22aμ, bν due to the lowest-order classical gluon fields in the McLerran-Venugopalan model of heavy ion collisions and observe an interesting transition from the classical fields to the quasi-particle picture at later times. Motivated by this observation, we push the formalism to higher orders in the coupling and calculate the contribution to G22aμ, bν coming from the diagrams representing a single rescattering between two of the produced gluons. We assume that the two gluons go on mass shell both before and after the rescattering. The result of our calculation depends on which region of integration over the proper time of the rescattering τZ gives the correct correlation function at late proper time τ when the gluon distribution is measured. For (i) τZ ≫ 1/Qs and τ − τZ ≫ 1/Qs (with Qs the saturation scale) we obtain the same results as from the Boltzmann equation. For (ii) τ − τZ ≫ τZ ≫ 1/Qs we end up with a result very different from kinetic theory and consistent with a picture of “free-streaming” particles. Due to the approximations made, our calculation is too coarse to indicate whether the region (i) or (ii) is the correct one: to resolve this controversy, we shall present a detailed diagrammatic calculation of the rescattering correction in the φ4 theory in the second paper of this duplex.


Introduction
The ultimate goal of heavy-ion collisions is to produce and study quark-gluon plasma (QGP). QGP is believed to be the primordial matter in our early universe after the Big bang and before the formation of nucleons. Heavy-ion collision experiments at RHIC and LHC provide us a golden opportunity to study such a new form of matter. The detectors only measure the properties of the system at a very late proper time τ ∼ 10 15 fm/c after the collision. (Here τ = √ t 2 − z 2 , with t being time and z labeling the collision axis.) QGP is believed to exist only in the first 10-20 fm/c after the collision: for the majority of the JHEP03(2018)158 remaining time the system reduces to a multitude of hadrons free-streaming towards the detector. The entire evolution history of bulk matter in the first 10-20 fm/c can be studied only by comparing theoretical calculations to the experimental data at roughly 10 15 fm/c. However, a consistent first-principles QCD formalism allowing to study the system from the very beginning of the collision to a moderate late time is still missing.
Hydrodynamic models could give a good description of the collective behavior seen in the experimental data (see [1] for a recent review). Hydrodynamics can be taken as an effective theory of the underlying quantum field theory near (local) thermal equilibrium [2]. Since bulk matter in heavy-ion collisions is far from thermal equilibrium at the very early stage, hydrodynamics breaks down at those early times. In practice, hydrodynamic models are switched on at some initial time τ 0 with the initial condition provided by other theoretical studies. In order to study the early stage of the collision one needs to employ the underlying field theory.
At the very early stage of a collision, a large number of saturated gluons are believed to be freed from the two nuclear wave functions (see [3,4] for a comprehensive review of the saturation/Color Glass Condensate (CGC) physics). In this case the classical Yang-Mills theory applies. It has been extensively studied in [5][6][7][8][9][10][11][12][13][14][15][16][17]. However, this approach can only give a highly anisotropized energy-momentum tensor with the ratio of the longitudinal to transverse pressures P L /P T approaching zero at later times [11-13, 15, 16]. Early pressure isotropization has been observed if certain types of vacuum quantum fluctuations are included in the classical field simulation [17,18]. In this case one has to deal with the dependence of the medium energy-momentum tensor on the lattice spacing [17,19]. This is a consequence of the non-renormalizability of the classical field approach with vacuum quantum fluctuations [20,21].
The Boltzmann equation has also been broadly used in heavy-ion collisions. It can be derived from two-point Green functions in quantum theory using the so-called quasi-particle approximation near thermal equilibrium [22][23][24][25][26]. The transition from classical fields to quasi-particles is expected to occur at τ ∼ 1/Q s with Q s the saturation momentum of the colliding nuclei [27]. Then, a parametric estimate using the quasi-particle picture gives a bottom-up scenario for the system to establish thermal equilibrium [27]. This picture has recently been confirmed by numerical solutions of the Boltzmann equation [28]. One of the intriguing questions about the Boltzmann equation is when it starts to apply to heavy-ion collisions since the derivation of this equation from quantum field theory has mostly been done for the systems near thermal equilibrium. The conventional understanding is that when the gluon density f is less than 1/g 2 (with g the QCD coupling), both the Boltzmann equation and classical field approximation apply [29]. However, this argument is based on the so-called quasi-particle approximation. It is of great interest to understand whether and how such a transition occurs in the collision process.
The Schwinger-Keldysh formalism or the close-time path formalism was first invented by Schwinger [30] and Keldysh [31]. It gives a unified description of equilibrium and nonequilibrium systems in quantum field theory [23]. This formalism has been used to study thermal equilibrium systems in thermal field theory [32,33]. It has also been used to study non-equilibrium systems by resumming two-particle-irreducible (2PI) or n-particle-JHEP03(2018)158 irreducible (nPI) diagrams [24,34,35]. The interested reader is referred to [24,36] for a comprehensive review of the nPI effective theories. A 2PI non-Abelian gauge theory would be of great interest to heavy-ion physics. However, the truncated 2PI effective action leads to gauge-dependent results for most observables [37]. In high-energy nuclear physics, the Schwinger-Keldysh formalism has been employed to resum leading ln 1 x terms with x the energy fraction into the color charge density functionals describing the colliding nuclei [38][39][40]. However, contributions beyond the leading ln 1 x have not been evaluated: such contributions could be important for the evolution of the system at late times.
The main purpose of this paper is to adapt the Schwinger-Keldysh formalism to study heavy-ion collisions in a perturbative approach. This approach is obviously gauge invariant. To show the advantages of this approach, we use the Schwinger-Keldysh formalism to try to derive the kinetic theory description of heavy ion collisions by studying a single 2 → 2 rescattering of the gluons produced in the collision in the saturation framework. To simplify the otherwise almost intractable calculation, we have to put the gluons on mass shell "by hand": this reflects the general physical expectation that the produced gluons would live long enough to go on mass shell before rescattering. This is a standard assumption in kinetic theory and should help us in deriving this theory. Making the onshell approximation requires us to assume that the proper time of the rescattering, τ Z , is long. In addition, in order to obtain kinetic theory, we need to assume that after the 2 → 2 rescattering the gluons again live long enough to go on mass shell, such that the observer at the proper time τ measures physical on-shell gluons. This implies that the time τ − τ Z should also be long, but, with the precision of our calculation, does not tell us which one of the two times, τ Z or τ − τ Z is (much) longer than the other or whether the two times are comparable. Unfortunately this turns out to be a crucial question in our calculation: if we assume that (i) τ Z 1/Q s and τ − τ Z 1/Q s (with both times comparable) and integrate over τ Z in this region, we reproduce the kinetic theory; however, imposing a particular ordering between the two time scales, namely using (ii) τ − τ Z τ Z 1/Q s in our τ Z integral gives us the so-called "free-streaming" situation, in which all gluons simply fly away without interacting. (The physical interpretation of the condition (ii) is not clear at present: it may represent some sort of a curious self-similarity pattern in the typical interaction times, where (τ − τ Z )/τ Z ∼ τ Z /(1/Q s ).) We conclude that a more detailed calculation without the on-shell assumptions is needed to determine whether the region (i) or (ii) (or neither) give us the correct late-τ asymptotics of the diagrams at hand. Such calculation, in the framework of the scalar ϕ 4 theory is carried out in [41]. This paper is organized as follows. We first give a brief review of this formalism and derive the Feynman rules for perturbative calculations in section 2. In section 3 we recalculate the gluon two-point function by using the lowest-order classical gluon fields of the McLerran-Venugopalan (MV) model [5][6][7] in the light-cone gauge. Based on this calculation, we show explicitly how quasi-particles emerge from classical fields. In section 4 we apply the Schwinger-Keldysh formalism to study the contribution from a 2 → 2 rescattering of these quasi-particles to the two-point Green function. That is, we study the rescattering of two particles produced by the classical gluon fields, assuming, as mentioned above, that the particles go on mass-shell both before and after the collision. The exact diagrammatic calculation at this order of the coupling in QCD is prohibitively complicated: instead we only keep the diagrams with the 2 → 2 rescattering and evaluate them using the on-shell approximation, by assuming that particles travel for long enough time to go on mass shell before and after the rescattering. This approximation is done in an effort to obtain the Boltzmann equation description of heavy ion collisions. The result of this calculation at late proper time τ (when the gluon is "measured") appears to depend on the region of integration over the rescattering proper time τ Z . For (i) τ Z 1/Q s and τ − τ Z 1/Q s our diagrammatic approach leads to the same answer as that obtained by solving the Boltzmann equation. However, as we show in section 5, for (ii) τ − τ Z τ Z 1/Q s the result is consistent with free-streaming gluons in the final state, and is very different from the solution of the Boltzmann equation. Further discussion of the physics behind the differences of cases (i) and (ii) is presented in section 6. The resolution of the question of whether the assumption (i) or assumption (ii) is correct is done in the second paper [41] of this duplex.

The Schwinger-Keldysh formalism for heavy ion collisions
In this section, we shall give a detailed description of the formalism used in our calculation. We adapt the Schwinger-Keldysh formalism [30,31] to describe the collision of two particles composed of a finite number of constituents. Following [5,6,8,42], the two colliding nuclei are taken to consist respectively of A 1 and A 2 constituent quarks at t = −∞, each valence quark representing a nucleon. Classical gluon fields resum the parameters α 2 s A to all orders [43]: in the actual calculations below we assume these parameters to be small, which would allow us to expand in them perturbatively.

The Schwinger-Keldysh formalism in perturbation theory
We formulate our problem in terms of the density matrix ρ, which can be written in terms of the wave functions of the two colliding particles Ψ 1 and Ψ 2 before the collision In the Schrödinger picture, the expectation value of any operator O is given by

JHEP03(2018)158
In order to perform perturbative calculations, we shall use the interaction picture by separating H into a free part H 0 and an interaction part H I . Let us denote the operator in the interaction picture by and the time evolution operator by where L I is the interaction Lagrangian corresponding to H I . From (2.2), it is easy to show that It is convenient to define the time ordering T C along the Schwinger-Keldysh contour C = C + C − . As illustrated in figure 1, the contour C runs from −∞ to +∞ and back to −∞. On C, the time ordering T C can be defined in the same way as the normal time ordering by replacing the θ function by [44] Accordingly, one can write Given any free field Φ, one can define the free propagator where Φ has been decomposed into positive-and negative-frequency parts, Φ p and Φ n .
Here F = 1 for fermions and F = 0 for bosons. Using this definition one can easily generalize Wick's theorem (see, for example, [45]) to the case of contour C by induction, that is, where the normal ordering operator N [· · · ] puts the negative-frequency parts to the left of all the positive-frequency parts in the product and each contraction of two fields gives JHEP03(2018)158 . The perturbative series can be generated by using series expansion of the exponential function in (2.5) With Wick's theorem in (2.9) and the propagator in (2.8), the above equation allows one to calculate O(t) perturbatively. Fields which are not contracted with other fields are to be contracted with either Ψ 1 , Ψ 2 | or |Ψ 1 , Ψ 2 .

QCD on the Schwinger-Keldysh contour
With the gauge fixing term, the QCD Lagrangian in n · A = A + = 0 light-cone gauge takes the form L can be separated into the free part and the interaction part In perturbative calculations, it is convenient to write the time integration of L I over C in (2.5) as where the field Φ represents any field in L I and the subscripts ± stand for the field Φ on C + and C − respectively. In the same notation as [20], we shall use the retarded/advanced basis in terms of the following fields with Φ ± the fields respectively on C + and C − contour. That is, in order to avoid dealing with the integration over C, one can double the number of fields instead. Accordingly, the propagator G(x 1 , x 2 ) can be taken as a 2 × 2 matrix in the space of the 1, 2 labels. In momentum space, for a scalar field with mass m we have

JHEP03(2018)158
For the quark field the propagator is (i, j are color indices) (2.18) and, for the gluon field in the limit ξ → 0, Perturbative calculations in QCD can be carried out using the interaction Lagrangian [40] with A µ ≡ A µ 2 and η µ ≡ A µ 1 . In summary, perturbative calculations of any operator can be carried out in momentum space by the following steps: 1. Draw all the Feynman diagrams at a certain order in g using the QCD vertices in (2.14).
2. Assign "1"s and "2"s to the fields at each vertex. All the allowed assignments have either one or three "1" fields at each vertex (see (2.20)). Keep in mind that (a) the contraction of any two "1" fields is always zero; and (b) the incoming states in the wave function of the colliding particles are only contracted with "2" fields. Therefore, each external line is assigned an index "2". The contraction results in spinors for external quarks and polarization vectors for external gluons in agreement with the conventional perturbative QCD.
3. Each contraction of any two fields gives the free propagator as one of the matrix elements of either (2.18) or (2.19). Here, the assignment of "1" and "2" gives the indices of the matrix element.
4. Each vertex is given by the corresponding one in the conventional perturbative QCD (see, say, [4]) with an overall prefactor 1/2 n 1 −1 with n 1 the number of "1" fields in this vertex.

5.
There is a conservation of 4-momentum at each vertex.
6. Integrate over each undetermined loop momentum.
7. figure out the overall symmetric factor of each diagram with a given assignment of "1"s and "2"s.
The above Feynman rules from steps 1 and 2 can be also obtained directly by using the Lagrangian with the doubled fields in (2.20).

JHEP03(2018)158
2.3 Modeling the nuclear wave function at t = −∞ To describe heavy ion collisions we need to augment the above Feynman rules by a specific definition of the density matrix. In this subsection, we take the same nuclear wave function at t = −∞ as those in refs. [5,7,8,42]. Big nuclei are taken to be composed of valence quarks at t = −∞. These quarks are confined in nucleons, which are homogeneously distributed inside the nuclei with a radius R. We shall study the collision of two big nuclei in the center-of-mass frame. Partons from nucleus 1 and 2 respectively have a large "+" and "−" , that is, these partons are approximately moving along their respective light-cones. The two nuclear wave functions are products of the wave functions of nucleons, which, in turn, are products of the valence quark wave functions.
The density matrix is We take the contribution to the density matrix coming from the "+" moving nucleus A 1 and write with the valence quark states |p i , p + i , m i . Here m i , n i are the quark color indices: summation is assumed over repeated indices. Define the Wigner distribution of a valence quark from nucleon i in nucleus A 1 by (cf. [46,47]) Substituting eq. (2.23) back into eq. (2.22) we obtain where P i = (p i + p i )/2. We have also approximated p + i ≈ p + i ≈ P + i since all the valence quarks in a relativistic nucleus have approximately the same light-cone momenta.
The averaging of an operatorÔ in the state |A 1 gives
In the standard MV model for a large unpolarized nucleus one usually neglects the transverse momenta P i of the valence quarks in the nucleons and assumes that the longitudinal momentum of the nucleus is evenly distributed among the nucleons. The corresponding quasi-classical Wigner function in the MV model is [46] with the nucleon number density ρ(b − , b) normalized such that and P + the light-cone momentum of the entire nucleus. Substituting eq. (2.27) into eq. (2.25) we arrive at where ρ 1 is the nucleon number density in nucleus A 1 . We have also suppressed the momenta in the argument of O in eq. (2.29): it is understood that P i = 0 and P + i = P + /A 1 for all the nucleons (or valence quarks) in the nucleus A 1 .
Since for a large nucleus in the MV model |A 1 , A 2 = |A 1 ⊗ |A 2 we conclude that the average over the initial states of a given operator O, which may represent a Feynman diagram, is are the positions of valence quarks in the nucleus A 2 while ρ 2 is the nucleon number density in that nucleus.
We see that the averaging in the nuclear wave functions in the MV model amounts only to averaging over positions and colors of the valence quarks in the two colliding nuclei [5][6][7][8].
In the following calculations, which are leading-order in A 1 and A 2 since they involve only one nucleon out of each nucleus, for simplicity we will put We will assume that the nuclei are identical, A 1 = A 2 = A, and have the same radii. Here R is the transverse radius of the nuclei and S ⊥ = πR 2 is the transverse cross-sectional area. propagators, crossed by orange dashed lines, separate each diagram into two. In A + = 0 gauge, there are three diagrams for the classical field A aµ . Each diagram in this figure corresponds to that in the product of two classical fields. In each diagram the quark on the top has a large P + while the one at the bottom has a large P − .

Classical fields and quasi-particles
In this section we calculate G aµ, bν   22 in the Wigner representation ) in the MV model power counting developed in [43]. In thermal field theory, the free correlation function is G 22 (X, p) = 2π(n B + 1/2) δ(p 2 ) with n B the Bose-Einstein distribution [44]. In systems near thermal equilibrium, one may neglect the dissipation near the quasi-particle peak in the spectral function and take G 22 (X, p) = 2π(f + 1/2) δ(p 2 ) with f the distribution function in order to derive the Boltzmann equation [23,24,26]. In this section, we shall study how the (quasi-)particle picture with p 2 = 0 emerges from the classical fields.

The classical field approximation at
In this subsection, we calculate G aµ,bν ). The lowest-order classical gluon field in covariant gauge was found before in [10,15]. G aµ,bν 22 is a gauge-dependent quantity. We shall show that it takes a much simpler form in A + = 0 gauge, which has a more transparent physical interpretation.

JHEP03(2018)158
We need only to evaluate the 9 diagrams 1 as shown in figure 2. In each diagram in this figure, the quarks are put on mass shell by each S (0) 22 propagator. As a result, each diagram corresponds to that in the product of two classical fields, in accordance with the discussion in appendix A. By including all the diagrams with possible crossing of internal gluon lines, we get (for the two identical nuclei described in section 2.3) where the trace, as defined in (2.30), puts l 2 = −l 1 . The classical field is Since we are interested in the mid-rapidity region, we only need the pole at k 2 = 0 and we can neglect the poles at k ± = 0. In this case, one can write By keeping only the logarithmically enhanced terms after integrating out l, that is terms with ln(k 1T /Λ) = ln(k 2T /Λ), we have and Λ is the infrared cutoff. This is much simpler than that in covariant gauge and we have checked that it gives exactly the same energy-momentum tensor as calculated in covariant gauge in [15] (see also [48,49]). 1 Here, we discard terms proportional to δ(x ± 1,2 ) in G aµ,bν

22
(x1, x2). Otherwise, there will be more diagrams. For example, one can not neglect the diagrams with the outgoing gluons attached to the quark at the bottom even in A + = 0 gauge when calculating the correlation function on the light cone.

From classical fields to quasi-particles
We write the retarded Green function in the following way where we dropped the i 's in all k + + i and replaced them by θ(x − ) in the prefactor. (The inverse Fourier transform would reinstate these i 's due to θ(x − ).) G aµ,bν 22 (x 1 , x 2 ) can be expressed in the following form Inserting the above expression into the first line of (3.1) and integrating out x gives To see the transition from the classical fields to quasi-particles, we need to study the late-time τ asymptotics of the correlation function. At large X + = τ e η / √ 2 and X − = τ e −η / √ 2 one is allowed to make the following approximations We get

JHEP03(2018)158
The two δ-functions give us two equations, which have two solutions are given by (3.14) Taking into account the above two solutions in (3.12) leads to At large X + , the predominant region of the above expression locates near p 2 0. By neglecting terms ∝ tan c X p 2 we have Our result can be further simplified by taking The above equation holds because the support of the left-hand side is limited to x = 0 as τ → ∞ and

JHEP03(2018)158
As a result, we have Our result in (3.21), while obtained in the classical field approximation, has a physical interpretation in terms of particles. We have taken the longitudinal size of the two nuclei to be zero in (2.31). As a result, they collide at t = 0 = z. After that, each produced gluon travels at the speed of light. Along the z-direction, its location From (3.22), one can easily see that the longitudinal pressure is zero at mid-rapidity due to δ(y − η), which indicates that the classical gluons simply free-stream away from the collision. This is what was earlier observed in [15]. Numerical simulations have shown that including all the other classical diagrams will not change the fact that the longitudinal pressure approaches zero much faster than the transverse pressure at late times [11][12][13][14]16], such that free-steaming is a property of the all-order classical gluon production in CGC, and not just of the lowest-order calculation carried out here.

Rescattering and the Boltzmann equation
In this section we will use G aµ,bν    Here, we only include the twoparticle irreducible (2PI) gluon self-energies in each diagram. Gain term: Loss term:  As illustrated in figure 4, we can group the subset of diagrams into a gain term and a loss term in the kinetic theory notation, with the circles denoting 2PI self-energies Π ij . In the gain term two classical G aµ,bµ In terms of the averaged self-energies the gain term takes the form To evaluate this expression we write the retarded Green function G 21 (x) in the form of (3.8), while the advanced and cut Green functions are and (4.5)

JHEP03(2018)158
Integrating over z, p , x we arrive at To reproduce kinetic theory one has to assume that gluons go on mass shell between interactions. This means the time between rescatterings is long enough for the gluons to go on mass shell. Therefore, we need to assume that X + − Z + and X − − Z − are very large in eq. (4.6). This approximation is different from simply assuming that X + and X − are large, as was done in eqs. (3.11), since the integrals over Z + and Z − in eq. (4.6) are not restricted to the regions far away from X + and X − respectively. Thus we simply assume that the large-X + − Z + and X − − Z − region dominates in the integral. This assumption is needed to obtain kinetic theory from our formalism, but cannot be easily justified otherwise for the collision at hand.
When assuming that a dimensionful quantity is large one has to compare it to another dimensionful quantity. Unfortunately this is hard in our case, since almost everything else is integrated out. We simply state that X + − Z + and X − − Z − are the largest distance scales in the problem, with the possible exception of Z + and Z − which may be comparable. Note that in deriving the classical correlator (3.21) we have assumed that X + and X − are large (see (3.11)): in the problem at hand, X + and X − from eq. (3.22) become Z + and Z − since we will be using the classical correlators to calculateΠ ij . Therefore, our Z + and Z − have already been assumed to be very large.
Finally, a question remains whether to send X + − Z + and X − − Z − to +∞ or to −∞ when taking them large: from the curly brackets in eq. (4.6) we see that only the X + − Z + → +∞ and X − − Z − → +∞ limits give a non-zero result. Applying those limits to eq. (4.6) with the help of eqs. (3.11) and integrating out k + , k + and Z afterwards while assuming thatΠ ij (Z + , Z − , Z, P ) =Π ij (Z + , Z − , P ) due to the slowly changing transverse JHEP03(2018)158 profile of the large nucleus yields where In arriving at eq. (4.7) we put p 2 = 0 in the argument of the Sign-function: this approximation will be justified shortly. Lower limits of the Z + and Z − integrals were set to zero in eq. (4.7) due to the classical correlator (3.22) which we will use to calculatē Π ij : the correlator ensures that no gluons are produced before the heavy ion collision at (t, z) = (0, 0). It is important to point out that, even though we assumed that X + −Z + and X − −Z − are very large, we have set the upper limits of the Z + and Z − integrations in eq. (4.7) to X + and X − respectively. This is related to the fact that our calculation requires that X + − Z + and X − − Z − are large, but does not tell us whether they need to be larger than Z + and Z − . For instance, large X + − Z + may imply either of the following situations (ditto for X − − Z − ): As we will see below, the two integration regions give different results. As we mentioned in the Introduction, we will answer the question of whether regime (i) or (ii) is correct at large X + , X − by a more detailed calculation in our next paper [41]. By assuming that τ X−Z is sufficiently large once again and using (3.19), we obtain Now let us turn our attention to the loss term. We will make similar approximations while evaluating the loss term in figure 4. The exact starting form of the loss term is as follows from (3.21) into (4.1). Similar to the above we define Z = (z 1 + z 2 )/2, z = z 1 − z 2 and writē Integrating over z, p and x yields where again y ± = x ± − z ± along with k = 2p − k. We have also assumed that x ± and z ± are much smaller than X ± and Z ± and neglected x ± and z ± in the argument of f cl .
Assuming that X + − Z + → +∞ and X − − Z − → +∞ andΠ ij (Z + , Z − , Z, P ) = Π ij (Z + , Z − , P ) we integrate over y + , y − , Z, k + and k + obtaininḡ where again we have put p 2 = 0 in the argument of the Sign-function along with the argument of f cl . Finally, invoking the late-time argument again we apply eq. (3.19) to eq. (4.15). This givesḠ . In each diagram the dashed orange line cuts through all the 2 − 2 propagators, separating the diagram into a product of an amplitude and a complex conjugate amplitude.

Gluon self-energies
We first evaluate the gluon self energy Π aµ,bν 11 from the diagrams in figure 5. Without loss of generality (for the late-time approximation at hand), we assume that the gluon propagators in these diagrams take the following form In our calculation, we choose to label by p 1 , p 2 , p 3 the momenta of the three 2 − 2 propagators in each diagram. These will be our integration variables. They satisfy p = p 1 + p 2 + p 3 . Each propagator has a positive-and negative-frequency part. We shall take the external momentum p + to be positive, and, in view of the above calculation of the gain and loss term, on mass shell, p 2 = 0. In each diagram, while evaluating loop integrals, there should be only two lines out of p 1 , p 2 , p 3 -carrying 2-2 lines with the positive-frequency parts of the propagators, while the remaining third line would come in with the negativefrequency part. It is clear that the diagrams in figure 5 reduce to the gg → gg scattering amplitude squared. Except for the first diagram in figure 5, different choices of positiveand negative-frequency parts for the 2-2 lines give us the products of s−, t− and u− channel amplitudes and their conjugates. Since p 1 , p 2 and p 3 are dummy variables to be integrated out, we redefine p 1 as the negative-frequency momentum and replace p 1 → −p 1 such that the new p 1 would have a positive frequency. Then, by collecting all terms obtained in this way, we arrive at the following result where and for brevity we have denoted with ω p = | p| and the Mandelstam variables defined by Next, let us calculate −i[Π 21 −Π 12 ], which are given by 2 times the real parts of the diagrams in figure 6. As indicated by the dashed lines in this figure, there are only two 2 − 2 propagators in each diagram. Compared to each corresponding diagram forΠ 11 in figure 5, the diagrams in figure 6 have a retarded (or advanced) propagator instead of the third 2 − 2 propagator. Then, by subtracting outΠ 12 fromΠ 21 one converts the retarded (advanced) propagator into a on-mass shell δ-function with different signs for its positiveand negative-frequency parts. After this, using the same trick as that for Π 11 with the positive and negative energy parts of the propagators, we get

Comparison with kinetic theory
In this subsection we evaluateḠ 22 (X, p) using the self-energies calculated in the previous subsection. Since both the gain (4.10) and loss (4.16) terms are proportional to δ(p 2 ), quasi-particle picture applies; therefore, as a comparison we also calculate the distribution function at O(α 2 s ) by performing a perturbative solution of the Boltzmann equation.

Results from the above approximation
Inserting (4.18) and (4.22) into (4.10) and (4.16) gives Both terms give a boost-invariant (rapidity-independent) particle distribution. This is more transparent if one uses the following variables Since our approximation should break down at early times, we require τ Z > τ 0 with τ 0 some initial time. Since f cl (X, p) ∝ δ(η − y), it is convenient to take Let us start with the gain term. Using eq. (3.22) we writē with p 23⊥ ≡ p 2⊥ + p 3⊥ . By integrating out y 1 and η Z we havē and the sum over {y 1 , η Z } goes over the following two values for each variable

JHEP03(2018)158
The above two solutions for {y 1 , η Z } respectively give Now let us integrate out τ Z . The δ-function in eq. (4.28) gives (4.32) and the Jacobian From the above equations, we finally obtain where η − η Z should have the same sign as y − η Z and y − η, and {y 1 , η Z } only assume the values in (4.30) which ensure that τ Z is positive.
For the loss term (4.24) evaluation appears to be more complicated in general. The difficulty is in the (X + Z)/2 in the argument of one of the f cl in (4.24). It appears that to obtain something similar to kinetic theory one has to replace in the argument of f cl in (4.24). This is an ad hoc assumption, particularly for the plus and minus components X ± , Z ± , which does not follow from the integration regions (i) and (ii) considered above. In fact, it can never be realized in the case of region (ii). Within the region (i) one could imagine a situation where and the replacement (4.35) may be justified. Such a condition is a further refinement of the integration region (i) and was not needed for the gain term. Below we will assume that the ordering (4.36) applies and make the substitution (4.35).

JHEP03(2018)158
Working in this approximation one can easily integrate out η Z and y 1 in eq. (4.24) and get At the end, we havē with {y 2 , y 3 } summed over two values each given by the solutions to the equations resulting from the last two δ functions in (4.37).
We are now going to compare the results (4.34) and (4.39) of this subsection for the gain and loss terms with the predictions of kinetic theory.

Results from the Boltzmann equation
In the boost-invariant and dilute (f 1) system, the Boltzmann equation (B.14) reduces to [50] (4.40) While the system we consider is boost-invariant (rapidity independent), we will concentrate on central rapidity, η = 0, throughout this subsection. At η = 0 one has τ = t. In order to see the connection to the calculation in our formalism, we write

JHEP03(2018)158
Knowing the distribution function f one can calculate the energy-momentum tensor using (4.47) Clearly, the initial conditions (4.42), or, equivalently, the classical gluon correlator (3.22) give T µν ∼ 1/τ for all the non-zero components of the energy-momentum tensor, along with the longitudinal pressure P L = T 33 (η = 0) = 0: this behavior corresponds to free streaming of gluons. Adding the correction f (1) from eq. (4.46) we obtain The exact values of the coefficients A (0) , A (1) and B (1) can be found by explicit integration: their exact values are not important to us, as long as A (1) and/or B (1) are not zero. The ∼ 1/τ and ∼ (1/τ ) ln(τ /τ 0 ) terms that A (1) and B (1) multiply constitute a deviation from the ∼ 1/τ free-streaming behavior of the classical gluon fields. We conclude that kinetic theory predicts a deviation from free streaming after including a single 2 → 2 rescattering correction to the classical gluon correlator. This prediction appears to agree with the results of the approximate calculations carried out above, after certain approximations were made. We will verify this prediction in [41].

Free streaming
Here we show how the above calculation can lead to different results depending on how the large-X + − Z + , X − − Z − condition is imposed. In particular, we demonstrate that the integration region (ii) from above does not lead to kinetic theory, but rather to free streaming of the produced gluons.

Free streaming after rescattering
In with a small parameter ξ 1 (but not too small, such that ξτ τ 0 still). For instance, we may have ξ = α λ s with some positive power λ > 0. Replacing τ → ξτ in the integration
It appears natural that in addition to the ordering (5.1) one also considers where ζ 1 is a large parameter. For τ 0 = 1/Q s one may have ζ = 1/α ω s with ω > 0 such that ζτ 0 = 1/(α ω s Q s ). While ζ is large, it should not be too large, such that ζ τ 0 τ still. eq. (5.2) is consistent with the condition (ii) from section 4.1 and also provides a way to impose the large-X + − Z + , X − − Z − condition.
Replacing the upper limit of the τ Z integration in eq. (4.28) by using τ → ζ τ 0 (along with the same replacement in the arguments of θ-functions, which we discard below since for small τ Z they are automatically satisfied) one gets We have employed τ z ζ τ 0 τ condition in simplifying eq. (5.3). (Strictly-speaking we have assumed a somewhat stronger condition Z ± ζ Z ± 0 X ± consistent with the original ordering (ii) to approximate τ X−Z ≈ τ and η X−Z ≈ η.) We see that in the integration region (ii), given by (5.2), the gain contribution to the correlation function still scales as 1/τ , but now it is multiplied by δ(y − η): the δ-function leads to zero longitudinal pressure, making thisḠ gain 22 (X, p) consistent with free streaming. The loss term is treated similarly: performing the τ → ζ τ 0 replacement in eq. (4.37) gives which is also consistent with free streaming. 3 We conclude that while the calculations of section 4.3.1 appear to be consistent with kinetic theory if the ordering (i) is imposed via (5.1), the same calculations are consistent with free streaming if the ordering (ii) is imposed with the help of (5.2). Therefore, since at this level of calculational precision we can not say whether the integration region (i) or (ii) is correct, we can not tell whether our calculation supports kinetic theory or the free-streaming scenario advocated in [15].

A general argument for free streaming
For completeness, let us briefly recap the free-streaming argument from [15], but now for the correlation functionḠ 22 considered in this work. (In [15] the argument was applied to the energy-momentum tensor of the medium produced in heavy ion collisions.) In a general case, involving all the possible multiple interactions and rescatterings, one can still write theḠ 22 correlation function as a sum of the three diagrams in the top row of figure 5, but now without requiring that Π ij include 2PI diagrams only. The circles now denote any (connected) diagram. The resulting expression is the same as in eq. (4.3): Since allΠ ij are Lorentz-invariant, they can only be functions of k 2 , k 2 and k · k . Thē Π ij 's can also be functions of p · k ∼ k − , p · k ∼ k + , p · k ∼ k − , and p · k ∼ k + with p and p the momenta of the nucleons in nuclei A 1 and A 2 . In [15] a "dictionary" was established by going through a number of examples. According to this "dictionary" the Fourier transform of the correlation function into coordinate space, converts (modulo some prefactors and an index shift of the Bessel function resulting from the transform) with η the space-time rapidity. According to eq. (5.7a), the leading late-τ contribution comes from putting k 2 = k · k = k 2 = 0 in theΠ ij 's, with corrections to it being suppressed by powers of 1/τ . Let us illustrate this using theΠ 11 term in eq. (5.5). Putting k 2 = k · k = k 2 = 0 in Π 11 and integrating that term in (5.5) over k − and k − yields In arriving at eq. (5.8) we have also neglected the k ± , k ± dependence inΠ 11 : below we will briefly outline how this can be reinstated. Employing

JHEP03(2018)158
we obtainḠ Hence the leading contribution to the correlator is ∼ 1/τ , corresponding to free streaming, as the energy-momentum tensor that one would obtain from the correlator (5.11) would scale as T µν ∼ 1/τ . If there exist terms proportional to, say, Bjorken hydrodynamics [51], which has T µν ∼ 1/τ 4/3 , they would be subleading compared to the free-streaming term of eq. (5.11).
Including the k ± , k ± -dependent terms in the above calculation would only add spacetime rapidity dependence in the correlation function owing to eq. (5.7b), without changing the conclusion (5.11) about the late-time asymptotics. Finally, the argument applies analogously to theΠ 12 andΠ 21 terms in eq. (5.5).
The remaining question is whether the leading asymptotics happens to have a zero coefficient, that is, what if In [15] it was shown that for the energy-momentum tensor this is not the case. The leading late-time contribution was shown to be proportional to the particle (gluon) production cross section. Specifically, the energy density was shown to be Hence the leading free-streaming term is non-zero as long as one can define the multiplicity of produced gluons dN/d 2 k T dη d 2 b ⊥ , that is as long as perturbation theory holds. 4

Conclusions and outlook
In this paper, we adapted the Schwinger-Keldysh formalism to study heavy-ion collisions in a perturbative QCD approach. We calculated the gluon two-point correlation function G aµ,bν 3 ) in the lowest-order classical approximation of the MV model. We found 4 The distribution of produced gluons, dN/d 2 kT dη d 2 b ⊥ , requires an IR cutoff once collinear-divergent corrections are included: such cutoff cancels in the integral of eq. (5.13) such that the energy density is independent of the cutoff (see [52]).

JHEP03(2018)158
that at large τ the (quasi-) particle picture emerges from the classical field calculation at this order in the sense that the resulting Green function is proportional to δ(p 2 ) and its Lorentz structure is determined by the polarization sum over the physical gluon polarizations, This Green function is also proportional to δ(y − η)/τ , which is a tell-tale sign of free streaming. We thus reproduced a well-known result that pure classical Yang-Mills dynamics in heavy ion collisions leads to free-streaming gluons [11-13, 15, 16]. Equipped with the Schwinger-Keldysh formalism as our calculational tool, we tried tackling the all-important problem in heavy ion collisions: how do these free-streaming gluons turn into QGP? An exact calculation of the first quantum correction to the classical gluon production in CGC appears to be prohibitively complicated in QCD, even if the gluon production is taken at the lowest non-trivial order. Instead we chose a subset of the diagrams at this order (order-α 8 s A 3 ) which correspond to the 2 → 2 rescattering of the classically produced gluons, simply because they are more likely to lead us to kinetic theory than other quantum correction diagrams. Each of these diagrams includes two subdiagrams of O(α 3 s A 2 3 ), described by G aµ,bν 22 in (6.1) each, along with the order-α 2 s Born-level cross section for the 2 → 2 scattering.
Even reducing the scope of the calculation to a subset of diagrams at the particular order in the coupling still did not make it tractable, and further simplifications turned out to be necessary. In order to obtain kinetic theory one needs to assume that the gluons travel for a long enough time to go on mass shell both before and after the rescattering. This consideration led us to making another ad hoc assumption that the time of the rescattering (τ Z ) and the time between the rescattering and the gluon "measurement" (τ −τ Z ) were both large, much larger than the gluon formation time of 1/Q s in the CGC. Note that τ is always assumed to be large in our calculation, since we are interested in the late-time dynamics: however, τ Z is an integration variable, and assumptions on τ Z and τ − τ Z imply restrictions on the τ Z integration range. Making the τ Z , τ − τ Z 1/Q s approximations allowed us to put the gluons on mass shell both before and after the collision, in complete analogy with the kinetic theory. However, it was still insufficient to obtain the Boltzmann equation.
To derive the Boltzmann equation we had to assume that while (i) τ Z 1/Q s and τ − τ Z 1/Q s , we do not impose any ordering between τ Z and τ − τ Z . This assumption, with a couple of others minor approximations, led us to the Boltzmann equation.
On the other hand, imposing a distinct ordering (ii) τ − τ Z τ Z 1/Q s between the two time scales in the τ z integral gave us free-streaming gluons with the correlation function similar to that in eq. (6.1), as it was also proportional to δ(y − η)/τ . We saw that after the rescattering the gluons still may assume a distribution similar to that for free-streaming particles. It seems surprising that after making several ad hoc assumptions designed to derive the kinetic theory we still could not rule out the possibility of the free-streaming gluons instead.
Our calculation above is an illustration of applying the Schwinger-Keldysh formalism to heavy ion collisions. Despite the simplifications brought about by the formalism, we needed to make a number of additional simplifying assumptions to complete the calculation. Our conclusion appears to be that we cannot tell whether a single rescattering of two classically produced CGC gluons brings about the first isotropization correction that one would expect for this process from the kinetic theory, as described in section 4.3.2. Our calculation neither fully confirms nor disproves the kinetic theory. We need to do a more careful diagram evaluation in order to see what the conclusion about the validity of kinetic theory in heavy ion collisions is. Unfortunately a complete diagram evaluation, free from all the above-listed assumptions, appears to be very complicated in QCD.
Luckily there is a simplification in the case of scalar field theories. In the companion paper [41] we perform a detailed calculation in the framework of the λϕ 4 theory to explicitly identify which approximation, (i) or (ii), is correct. The calculation in [41] leads to freestreaming scalar particles, thus indicating that the condition (ii) is correct, at least in the scalar theory, thus challenging the applicability of the kinetic theory to heavy ion physics.

JHEP03(2018)158
lines. Specifically, we shall show that in this limit the gluon two-point function G aµ,bν 22 reduces to a product of classical gluon fields, i.e., with the angle brackets denoting the averaging from eq. (2.30). As illustrated in figures 7(a) and 7(b), we shall prove that any diagram of order 1 3 ) nq that does not vanish in the eikonal approximation is separated into two disconnected pieces by the set of all S (one on each quark line) in the eikonal approximation. Each piece is connected by retarded Green function and hence contributes to the classical field.

A.1 The eikonal approximation of the quark lines
At high energies, the recoil of the valence quarks from radiating soft gluons is negligible. In this case one can make the so-called eikonal approximation to the quark lines [5,7,8,42]. In each diagram for our problem, the valence quarks of nucleus 1 receive some momentum transfer, l, from the scattering with other partons. l is typically much softer than P + 1 . This allows us to approximate the free quark propagator by . (A.2) Here we ignore the quark mass. Similarly, for a valence quark of nucleus 2 with momentum transfer l, one has Here, we have kept / l because it may not always give a suppressed contribution compared to / P 2 in the A + = 0 light-cone gauge.
In this paper we shall not consider the small-x quark production. That is, all the quark lines come from the two nuclear wave functions. In this case the eikonal approximation only involves replacing the quark propagator in (2.18) by eqs. (A.2) and (A.3). Below we shall show that substituting these two equations for each quark propagator S in the diagrams of order 1 g 2 (g 4 A

A.2 Diagrams in the classical field approximation
Let us focus on a generic diagram with n q valence quark lines. First of all, the diagram should be connected. Otherwise, it can be separated into the product of connected subdiagrams. Among these sub-diagrams, there must exist at least one diagram either without any gluon radiation or with one gluon radiation. Due to unitarity, connected diagrams without radiation cancel. And diagrams with one radiated gluons also vanish after the average over the initial distribution in (2.30) due to the color neutrality of the two-nuclei source. Therefore, it has to be connected.

JHEP03(2018)158
Second of all, the diagram should be of order g 4nq−2 A nq 3 in order to give a non-vanishing contribution in the classical limit. The factor A nq 3 ∝ A S ⊥ nq results from the average at the initial time in eq. (2.30). We shall prove that the diagram gives a non-vanishing contribution to G aµ,bν 22 in the eikonal approximation only if it has a S (0) 22 on each quark line. An example with n q = 4 is shown in figure 7. We will show that the first two diagrams, each of which is separated into two sub-diagrams connected by S In order to prove the above statements, we need first to count the number of 2-2 propagators in the diagram. Let us assume that there are n cl 4 four-gluon vertices with one "1" field, n qu 4 4-gluon vertices with three "1" fields, n cl 3 three-parton vertices with one "1" field and n qu 3 three-parton vertices with three "1" fields. Then, the diagram is of order g n D with with n qu = n qu 3 + n qu 4 the number of (quantum) vertices with three "1" fields. The second ingredient of our proof is that the diagram gives a non-vanishing contribution in the eikonal approximation only if there is at least one S (0) 22 on each quark line. Otherwise, its contribution will be canceled by other diagrams. Let us only single out one quark line without any S (0) 22 propagator. Assume that there are n gluon lines connected to it. First, we take n = 2. Let us include the diagram with the two gluon lines connected in the opposite order (while the rest part of the diagram is kept the same). The two diagrams cancel with each other due to the following cancellation with l ∓ 1 respectively corresponding to the diagrams with the quark being from nucleus 1 or 2. Here, we have used the fact that l ± 1 + l ± 2 ≈ 0 in the eikonal limit.

JHEP03(2018)158
1. The diagrams for the two-point gluon correlator G µν 22 should be of order g n D with n D ≥ 4n q − 2.
Indeed, eq. (A.6) with n ext 2 = 2 gives n D = 2n 22 + 2n q − 2 + 4n qu ≥ 4n q − 2 + 4n qu ≥ 4n q − 2, (A.12) where we have used the fact that n 22 ≥ n q since, for the diagram not to cancel, each valence quark line should contain a 2-2 propagator due to the proof above. The lowest possible value of n D corresponds to the classical dynamics. It is reached if n 22 = n q and n qu = 0. The latter condition means no vertices with three "1" fields in the diagrams for the classical correlator. This is consistent with the conclusion in the functional approach [29]. Each quark line has one S 22 , which separates the diagram into two sub-diagrams connected by the cut quark propagators. We conclude that classical G µν 22 must be a product of two classical gluon fields.

B The Boltzmann equation for gluons
In this appendix we review the standard derivation of Boltzmann equation for gluons. Let us define Inserting (4.18) and (4.22) into the above equation and ignoring 1 2 associated with f in (B.10) and the above equation gives the Boltzmann equation in the classical limit (see [53] for another derivation). If one keeps 1 2 's [29], one has (B.14) Just like in the λϕ 4 theory [21], the last term on the right-hand side of (B.14) leads to a UV divergence. This qualitatively helps understand the origin for the lattice spacing dependence observed in [17] although a quantitative analysis requires the calculation in lattice QCD (see a discussion in QED [54,55]).
In order to cancel the above UV divergence, one needs to include the contributions from the diagrams in figures 8 and 9. The diagrams in figure 8 give with |M | 2 given in (4.19).
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.