Time-Dependent Observables in Heavy Ion Collisions 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 $G_{22}^{a\mu,b\nu}$ 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 $G_{22}^{a\mu,b\nu}$ 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 the ordering between the proper time of the rescattering $\tau_Z$ and the proper time $\tau$ when the gluon distribution is measured. For (i) $\tau_Z\gg 1/Q_s$ and $\tau-\tau_Z\gg 1/Q_s$ (with $Q_s$ the saturation scale) we obtain the same results as from the Boltzmann equation. For (ii) $\tau-\tau_Z\gg \tau_Z\gg 1/Q_s$ 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 ordering (i) or (ii) is the correct one: to resolve this controversy, we shall present a detailed diagrammatic calculation of the rescattering correction in the $\varphi^4$ theory in the second paper of this duplex.

1 Introduction The ultimate goal of heavy-ion collision 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 time τ ∼ 10 15 fm/c after the collision. QGP is believed to exist only in the first 10-20 fm/c after the collision: for the majority of the 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). 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 energymomentum tensor on the lattice spacing [17,19]. This is a consequence of the nonrenormalizability 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 quasiparticle 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 non-equilibrium 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-irreducible (nPI) diagrams [24,34,35]. The interested reader is referred to [24,36] for a comprehensive review for 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 highenergy 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. This paper is organized as follows. We first give a brief review of this formalism and derive the Feynman rules for perturbative calculations in Sec. 2. In Sec. 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 Sec. 4 we 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 that the particles go on mass-shell both before and after the collision. The result of this calculation appears to depend on the ordering between the rescattering proper time τ Z and the proper time τ when the gluon is measured. 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 Sec. 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 Sec. 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 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 Fig. 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 G (0) (x 1 , x 2 ). 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 (2.13) 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 .
For the quark field the propagator is (i, j are color indices) 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).

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 , k . Here k i , l 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]) 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 and the curly brackets in the argument imply dependence on all the momenta or coordinates, e.g., {P i } = P 1 , P 2 , . . . , P A 1 .
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 , and have the same radii. Here R is the transverse radius of the nuclei and S ⊥ = πR 2 is the transverse cross-sectional area.

Classical fields and quasi-particles
In this Section we calculate G aµ, bν   22 in the Wigner representation 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. 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.
We need only to evaluate the 9 diagrams 1 as shown in Fig. 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 Sec. 2.3) 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.
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 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]).

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 At large X + and X − one is allowed to make the following approximations We get 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 τ 1 and 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 − η). This is what has been 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].

Rescattering and the Boltzmann equation
In this Section we will use G aµ,bν  3 ). This subset of diagrams can be obtained by assigning "1"'s and "2"'s to each diagram in Fig. 3 and replacing two of its 2-2 propagators with the classical one in (3.21). That is, the two of the 2-2 propagators are replaced by the 9 diagrams in Fig. 2 with all the possible crossings of their internal gluon lines. We shall show that under a certain approximation these diagrams give a result identical to that obtained by solving the Boltzmann equation via perturbative expansion in the collision term. In this sense they give the contribution to G aµ,bν 22 (X, p) from rescattering between the produced gluons. However, under a different approximation these diagrams do not reduce to a solution of Boltzmann equation. Here, we only include the two-particle irreducible (2PI) gluon self-energies in each diagram.  For the simplicity of the color and Lorentz indices, we shall calculatē As illustrated in Fig. 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µ

22
are used in the calculation of the self-energies, Π's, while in the loss term only one classical G aµ,bµ 22 is used in Π's and the other classical correlator is placed on one of the external gluon propagators, as shown by the green ovals in Fig. 4.
In terms of the averaged self-energies the gain term takes the form To evaluate this expression we write the retarded Green function G and 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 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 calculateΠ 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 limits give different results. As we mentioned in the Introduction, we will answer the question of whether regime (i) or (ii) is correct 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 Fig. 4. The exact starting form of the loss term is as followsḠ whereḠ cl 22 is obtained by substituting the classical correlator G aµ,bν 22 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

Gluon self-energies
We first evaluate the gluon self energy Π aµ,bν 11 from the diagrams in Fig. 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 with g 22 (X, p) = 1/2 for the free one and g 22 (X, p) = f cl (X, p) for the classical one. For our problem, we need only to include vertices with only one "1" field. In each diagram there are three 2 − 2 propagators according to the counting rule in (A.6). 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 negative-frequency part. It is clear that the diagrams in Fig. 5 reduce to the gg → gg scattering amplitude squared. Except for the first diagram in Fig. 5, different choices of positive-and 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  retarded (advanced) propagator into a on-mass shell δ-function with different signs for its positive-and 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ē where (4.29) and the sum over {y 1 , η Z } goes over the following two values for each variable 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 orderings (i) and (ii) considered above. In fact, it can never be realized in the case of ordering (ii). Within ordering (i) one could imagine a situation where and the replacement (4.35) may be justified. Such a condition is a further refinement of the ordering (ii) and was not needed for the gain term. Below we will assume that the ordering (4.36) applies and make the substitution (4.35).
Working in this approximation one can easily integrate out η Z and y 1 in Eq. (4.24) and get with η Z = y. Now we are left with 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 The terms f (n) are to be found from solving the Boltzmann equation order-by-order in the coupling α s . The initial condition (the value of f before the collision term becomes important) is given by saturation dynamics (cf. Eq. (3.22)), which satisfies the "free" Boltzmann equation The higher orders in α s can be calculated by iteration with C (n+1) the constant of integration at O(α 2(n+1) s ) and t 0 some initial time when the Boltzmann dynamics starts to apply (e.g. t 0 ∼ 1/Q s ).
For comparison with the results of the previous Subsection we need only to evaluate terms proportional to O(α 2 s ). According to (4.44), we first evaluate To satisfy the initial conditions at time t = t 0 we put the integration constant to zero, C (1) = 0. Substituting Eq. (4.45) into Eq. (4.44) and integrating yields We obtain a linear combination of ∼ 1/t and ∼ δ(y) (1/t) ln(t/t 0 ) terms. This is exactly the same t and y-dependence as that in the previous Subsection for the gain and loss terms respectively (if we apply η = 0 to those results).
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 − is imposed. In particular, we demonstrate that ordering (ii) from above does not lead to kinetic theory, but rather to free streaming of the produced gluons.

Free streaming after rescattering
In the calculations of Sec 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 limits of Eq. (4.27) would still lead toḠ gain 22 (X, p) ∼ 1/τ , just like in Eq. (4.34). Similarly, replacing τ → ξτ in the upper integration limit and in the theta-functions of Eq. (4.37) one still obtainsḠ loss 22 (X, p) ∼ (1/τ ) ln(τ /τ 0 ) δ(y − η), just like in (4.39). While the prefactors may be modified by the τ → ξτ substitution, the τ and η dependence of the gain and loss terms would remain the same. Hence, Eq. (5.1) appears to provide a more proper way of imposing the condition (i) on the calculation in Sec. 4.3.1.
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 Sec. 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 after 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 with the ordering (ii), (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 Sec. 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 ordering (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 Fig. 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 nucleus 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 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 space-time 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 that at large τ the (quasi-) particle picture emerges from the classical field calculation at this order in the sense that  4 3 ), which are beyond the classical field approximation, and corresponds to a 2 → 2 rescattering of the classically produced gluons. Each of these diagram includes two subdiagrams of O(g 6 A 2 3 ), described by G aµ,bν 22 in (6.1) each. In our calculation the rescattering occurs at some space-time point Z µ while the gluon distribution is measured at another point X µ . We made the following approximations: 1. τ Z ≡ √ 2Z + Z − 1/Q s Under this assumption, each sub-diagram took the form in (6.1).
1/Q s Under this assumption, the gluons after the rescattering can be taken to be quasiclassical particles. That is, they travel along a classical trajectory X 3 − Z 3 = (X 0 − Z 0 )p 3 /p 0 with p µ the four-momentum of the gluons.
Under the above approximations, which are equivalent to case (i) listed above (or, for the loss term, by additionally imposing τ Z τ X−Z ), we find the rescattering correction consistent with a power series in α s solution of the Boltzmann equation.
However, one needs to make a more detailed calculation to justify our approximations above. In our approximations of Eqs. (4.27) and (4.37) we put the upper limit of τ Z integration to be τ in apparent violation of the assumption 2. We also do not know a priori whether our assumption 2 correctly represents the full diagrammatic calculation, since one may take another limit instead, τ τ Z , corresponding to case (ii) by the above counting. In this caseḠ gain 22 andḠ loss 22 are still respectively given by (4.27) and (4.37) with another upper limit for τ Z integration, as detailed in Sec. 5. If one takes the limit τ → ∞ in this case, both the gain and loss terms are proportional to δ(y − η)/τ . That is, after the rescattering the gluons assume a distribution similar to that for freestreaming particles. In the companion paper [41] we perform a detailed calculation in the framework of the λϕ 4 theory to explicitly identify which approximation is correct.
22 cut quark propagators. Each of the sub-diagrams contributes to the classical field. In contrast, diagram (c) does not have such a separation. However, it is canceled in the eikonal approximation. Here, each orange dashed line crosses a quark 2 − 2 propagator indicating that the quark is on mass-shell and each retarded Green function is associated with an arrow pointing in the increasing time direction.

A The classical field limit
In this Appendix we will verify the classical field limit of our formalism. The limit involves (a) making the eikonal approximation of the quark lines from the nuclear wave functions; and (b) keeping only the diagrams of order 1 3 ) nq with n q the number of the quark 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 Figs. 7(a) and 7(b), we shall prove that any diagram of order 1 g 2 (g 4 A 1 3 ) nq that does not vanish in the eikonal approximation is separated into two disconnected pieces by the set of all S

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 sub-diagrams. 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. Second of all, the diagram should be of order g 4nq−2 A nq 3 in order to give a nonvanishing contribution in the classical limit. The factor A on each quark line. An example with n q = 4 is shown in Fig. 7. We will show that the first two diagrams, each of which is separated into two sub-diagrams connected by S (0) 22 on each quark line, give non-vanishing contributions while the third one vanishes in the eikonal approximation. Including all the non-vanishing diagrams such as those in Figs. 7 (a) and (b) leads to the classical field approximation in (A.1).
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.
This exactly cancels the other terms in I n+1 . By induction, (A.8) is true for all n. Finally, by using the power counting (A.6) and the identity (A.8), we can make the following statements about the classical field limit in our formalism: 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 In order to cancel the above UV divergence, one needs to include the contributions from the diagrams in Figs. 8 and 9. The diagrams in Fig. 8 give with |M | 2 given in (4.19).