Time-Dependent Observables in Heavy Ion Collisions II: in Search of Pressure Isotropization in the $\varphi^4$ Theory

To understand the dynamics of thermalization in heavy ion collisions in the perturbative framework it is essential to first find corrections to the free-streaming classical gluon fields of the McLerran-Venugopalan model. The corrections that lead to deviations from free streaming (and that dominate at late proper time) would provide evidence for the onset of isotropization (and, possibly, thermalization) of the produced medium. To find such corrections we calculate the late-time two-point Green function and the energy-momentum tensor due to a single $2 \to 2$ scattering process involving two classical fields. To make the calculation tractable we employ the scalar $\varphi^4$ theory instead of QCD. We compare our exact diagrammatic results for these quantities to those in kinetic theory and find disagreement between the two. The disagreement is in the dependence on the proper time $\tau$ and, for the case of the two-point function, is also in the dependence on the space-time rapidity $\eta$: the exact diagrammatic calculation is, in fact, consistent with the free streaming scenario. Kinetic theory predicts a build-up of longitudinal pressure, which, however, is not observed in the exact calculation. We conclude that we find no evidence for the beginning of the transition from the free-streaming classical fields to the kinetic theory description of the produced matter after a single $2 \to 2$ rescattering.


Introduction
The problem of isotropization and thermalization of the medium produced in ultrarelativistic heavy ion collisions is arguably the central theoretical problem in the field since it addresses the fundamental question of whether and how quark-gluon plasma (QGP) is formed in these collisions. Despite a number of theoretical efforts, the solution of this problem still remains elusive. Thermalization appears to be easier to tackle at strong ('t Hooft) coupling in the framework of the anti-de Sitter/Conformal Field Theory (AdS/CFT) correspondence [1,2]: there it is possible to show that a collision of two shock waves results in the black hole formation in the AdS 5 bulk, corresponding to a thermal medium being formed at the boundary. This was demonstrated analytically (but indirectly) using the trapped surface analysis [3][4][5] and was directly observed in a numerical solution of Einstein equations [6,7]. The weakness of the approach based on AdS/CFT correspondence is that the duality is for N = 4 super-Yang-Mills theory and not for quantum chromodynamics (QCD). Nevertheless, a consensus exists in the community that thermalization of the medium produced in high-energy collisions in strongly-coupled field theories is very likely to take place and to happen on a very short time scale.
The efforts to tackle the isotropization and thermalization problems at weak coupling have not achieved such a universal consensus. The initial break-through in the theoretical understanding of thermalization in QCD at weak coupling was the so-called 'bottom-up thermalization' proposal [8]. In this scenario, the quark-gluon system that begins in an initial state due to saturated gluon fields created in nuclear collisions (dominated by the classical gluon fields of the McLerran-Venugopalan (MV) model [9][10][11][12][13][14][15]) progresses to a thermalized isotropic medium due to 2 → 2, 2 → 3 and 3 → 2 rescatterings. However, the qualitative arguments presented in [8] have never been verified by explicit diagrammatic calculations. Moreover, in [16][17][18] it was pointed out that the bottom-up thermalization scenario may be invalidated by the occurrence of plasma instabilities, which could be present due to the momentum-space anisotropy of the initial non-equilibrium gluonic medium. Numerical simulations of these instabilities appear to demonstrate that the growth of instabilities is stopped due to the non-Abelian nature of strong interactions [19,20], possibly reinstating the original 'bottom-up' scenario. Alternative approaches [21] apply classical Yang-Mills dynamics to confirm the parametric scaling of the observables predicted in the first stage of the 'bottom-up' thermalization [8], and possibly leading to eventual thermalization of the medium [22]: however, the resulting classical dynamics appears to be non-renormalizable [23,24].
Yet another weakly-coupled approach to thermalization is based on Boltzmann equation. The applicability of Boltzmann equation to description of the medium produced in late stages of heavy ion collisions was argued in [25,26] (with the Vlasov-Boltzmann equation used in the instability studies mentioned above). Boltzmann equation dynamics appears to lead to thermalization of the produced medium [27], confirming the bottom-up thermalization scenario. However, the existing literature lacks a side-by-side comparison of Boltzmann equation with the explicit diagram calculation for heavy ion collisions: such a comparison is needed to either validate the applicability of the Boltzmann equation to heavy ion collisions or to prove otherwise. Performing such a cross-check is the main goal of the present work. In [25] a correspondence was established between the Boltzmann equation containing only the order-f 3 part of the collision term and the classical gluon fields at late times (and, hence, the late-time limit of the Feynman diagrams describing the collision in the classical approximation). However, the correspondence was never checked for the Boltzmann equation with the order-f 2 part of the collision term and the diagrams describing some correction to the classical fields of the MV model. This will be performed below.
In more general terms, many perturbative thermalization scenarios assume that the classical gluon fields of the MV model become sub-leading at late proper times in the collisions, being superseded by fields produced by some other dynamics, for instance due to the Boltzmann equation. However, no explicit calculation of Feynman diagrams exists in the literature which starts with the actual collision of two large nuclei, identifies a particular diagrammatic correction to the classical gluon fields and shows explicitly that such a correction becomes dominant at late times. The absence of such calculations is probably attributable to their complexity. However, an approach like this would have been natural in the saturation/Color Glass Condensate (CGC) framework [28][29][30][31][32][33][34]. There (and elsewhere in perturbative calculations in field theory) one usually starts with the tree-level leading-order contribution, which is often classical. The leading-order contribution receives corrections due to quantum fluctuations, leading parts of which may be resummed using evolution equations. This program has been carried out for the case of deep inelastic scattering (DIS) at small Bjorken x, where the leading order contribution to unpolarized DIS structure functions is given by the Glauber-Mueller quasi-classical multiple rescatterings [35], and the quantum corrections resumming logarithms of 1/x are included via the Balitsky-Kovchegov (BK) [36][37][38][39] and Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) evolution equations [40][41][42][43].
For heavy ion collisions analyzed in the saturation framework the leading contribution to, say, the energy-momentum tensor of the produced medium is given by the classical gluon field of the MV model [9][10][11][12][13][14][15]. This is already a very difficult calculation, only possible to be fully done numerically due to the complexity of the analytic attempts [44,45] (see [46][47][48][49][50] for perturbative results valid for proton-proton and proton-nucleus collisions). The numerical calculations [12][13][14][15] indicate that the classical gluon fields lead to a free-streaming medium, characterized by zero longitudinal pressure P L = 0 and the energy density = 2 P T ∼ 1/τ at late proper times τ 1/Q s . (Here P T and P L are the transverse and longitudinal pressures at mid-rapidity, τ = √ t 2 − z 2 is the proper time, and Q s is the classical gluon saturation scale.) Quantum corrections to the classical energy-momentum tensor resumming leading logarithms of 1/x were addressed in [51,52], where it was argued that such corrections can be resummed using the JIMWLK evolution equation for the weight functionals of the color charge densities in the two nuclei. The resulting gluon fields are still obtained by solving the classical Yang-Mills equations, but now with the JIMWLK-modified distribution of color sources. Therefore, such small-x evolution corrections still lead to a classical free-streaming energy-momentum tensor and are not related to isotropization or thermalization of the medium.
The question of whether perturbative quantum corrections to the classical gluon fields which usher in isotropization and thermalization exist still remains open. Over a decade ago, one of the authors of this work tried looking for such corrections in [53] (see also [54]). Having failed to find them, he argued that such corrections do not exist as long as one can define a gluon production cross section: hence, the end state of any perturbative (weakly-coupled) dynamics in heavy ion collisions was argued to always be a free-streaming bunch of particles [53,54].
The arguments of [53,54] notwithstanding, potential candidates for the isotropizationinducing quantum corrections are the 2 → 2 rescatterings as resummed by Boltzmann equation in the framework of kinetic theory. In the previous part I of this paper duplex [55] we showed that if one starts with the Boltzmann distribution function f (0) for the classical gluon fields of the MV model, and inserts it into the order-f 2 part of the collision term of the Boltzmann equation, solving the latter for a corrected distribution f (1) , one indeed does obtain isotropization corrections to the P L = 0, ∼ 1/τ free-streaming behavior of the classical gluon medium. The remaining question is whether Boltzmann equation correctly represents the Feynman diagrams it purports to sum. In [55] we review the derivation of Boltzmann equation that exists in the literature, concentrating on the same case of a single 2 → 2 rescattering of two classical gluon fields. The conclusion reached in [55] is that the underlying Feynman diagrams including the 2 → 2 rescattering lead either to results consistent with Boltzmann equation prediction of to free-streaming depending on how the late-time limit is taken. Namely, denote by τ 0 the time of the 2 → 2 rescattering (assumed to be instantaneous in the derivation, which is valid at late times only when the gradient expansion becomes possible) and by τ the time in the argument of f (that is, the time when we measure the particle in question). For the Boltzmann equation to be valid, the particles (gluons) must approximately go on mass-shell both in the time after they are produced in a collision but before they rescatter and in the time after they rescatter but before they are detected. This means that time intervals τ 0 and τ − τ 0 should be sufficiently long. While it is clear that for τ 0 "sufficiently long" (on the average) means τ 0 1/Q s , since 1/Q s is the time it takes for the classical gluon fields to go on mass shell, it is less clear what "sufficiently long" means for τ − τ 0 . In [55] we consider two options, and show that the ordering (i) yields results consistent with the Boltzmann equation, while the ordering (ii) leads to free streaming and is not consistent with the Boltzmann equation. Unfortunately the calculation performed in [55] (along with the earlier arguments in favor of Boltzmann equation) was too coarse to tell us whether the ordering (i) or the ordering (ii) follows from the full Feynman diagram calculation. It is the goal of the present paper to resolve this ambiguity, at least in the framework of the ϕ 4 theory that we use for simplicity instead of QCD. Below we will calculate the Feynman diagrams contributing to the 2 → 2 rescattering of two classical fields using the Schwinger-Keldysh formalism. As we have already mentioned, we will use the scalar ϕ 4 theory coupled to an external current [56] for simplicity. Without going into detail of the scalar particle production (though one could think of the Higgs production via gluon fusion), we assume that two scalar particles were produced in the collision with their distribution given by the two-point correlation functions very similar to that for the classical gluon fields in the saturation/CGC physics. (These two particles do not have to be on mass shell.) The setup of the problem is presented in Sec. 2 below. The particles rescatter via the 2 → 2 process, which is simpler in the ϕ 4 theory than in QCD: for instance, this interaction is truly instantaneous in the ϕ 4 theory. The two-point coordinate-space correlation function G(x 1 , x 2 ) resulting from the rescattering is calculated in Sec. 3. A calculation of the mixed-representation Green function G(X, P ) (that is commonly used in derivations of Boltzmann equation) resulting from the same 2 → 2 rescattering process is presented in Sec. 4. (Sec. 3 also contains a calculation of the corresponding energy-momentum tensor.) Both Green function calculations lead to the result consistent with free streaming and hence with the case (ii). Therefore, we see no evidence supporting the use of Boltzmann equation in describing the (perturbative) dynamics of the medium produced in heavy ion collisions. Our conclusions are summarized in Sec. 5.
2 Isotropization problem for the ϕ 4 theory

Classical two-point correlation function
The two-point 22 correlation function due to the lowest-order classical gluon fields in the MV model was calculated in [55] using the A + = 0 light-cone gauge. The result is with the polarization vector µ λ (k) = 0, λ · k k + , λ and λ=± i λ * j Here we assume the collision of two identical large nuclei with atomic numbers A 1 = A 2 = A, each of them shaped as a longitudinally-oriented cylinder with a very large cross-sectional area S ⊥ . Underlined variables denote two-dimensional vectors in transverse plane, v = (v 1 , v 2 ), while the light-cone variables are v ± = (v 0 ± v 3 )/ √ 2 with x 3 = z the collision axis. The contributing diagrams for the correlator (2.1) are shown in Fig. 3 of [55] and are comprised of two sets of lowest-order gluon production diagrams (cf. [46][47][48]).
As described above, it would be very interesting and important to find the perturbative correction to the correlator (2.1) due to a 2 → 2 rescattering process involving two such Green functions. However, full calculation of the order-α 2 s correction to Eq. (2.1) involving two classical correlators (that is, a calculation of the order-α 8 s (A/S ⊥ ) 4 correlator) appears to be prohibitively complicated in QCD. Instead we will tackle a similar problem in massless ϕ 4 theory. To do so we first have to construct an analogue of the correlator (2.1) in the massless scalar theory. This is achieved by replacing the polarization sum in Eq. (2.1) by (−1) and writing the rest of the expression as with f (k T ) a function of the magnitude of the transverse momentum k T = |k| which falls off rather fast at large k T and is infrared (IR) finite due to saturation effects. The exact form of f (k T ) is not going to be important below. (This function is proportional to the k T spectrum dN/d 2 k T dy of the produced particles in the classical approximation [53].) The rapidity-independent correlation function (2.3) can not result from a collision of particles taken entirely in the scalar theory: instead, one can think of it as resulting from some gluon+gluon→scalar fusion process, with the two gluons coming from the classical fields of the two colliding nuclei, as schematically shown in Fig. 1, where the shaded circle represents an effective gluon+gluon→scalar vertex. Higgs production through gluon fusion (via a top-quark loop) is one example of such a process (though of course Higgs is massive, unlike the massless scalar considered here). Let us stress one more time that the exact origin of the correlation function (2.3) is not important here: what is important is that it carries the main features of the gluon correlation function (2.1).
Our notation for the correlation function (2.3) is shown in Fig. 2, where the Feynman diagrams contributing to the correlator are summarily shown by a green oval. k Figure 1. Scalar particle production as envisioned here, with the solid lines denoting colliding quarks in the two nuclei and the dashed line denoting the scalar particle. The shaded circle may represent a quark loop, as would be the case in the Higgs production. In coordinate space the Green function (2.3) is given by the Fourier transform Employing the integrals in the appendices of [53] we obtain Using the large-argument asymptotics of Bessel functions we conclude that where we have averaged the cosine squared over time. (We have also employed the fact that f (k T ) is regulated by the saturation scale Q s in the IR, such that the k T ≈ 0 region does not contribute significantly to the integral.) Importantly the classical correlation function scales as which is a tell-tale sign of free streaming. To understand this better, let us calculate the energy-momentum tensor corresponding to this classical dynamics. The correlation function (2.5) is independent of the space-time rapidities η 1 = (1/2) ln(x + 1 /x − 1 ) and . It is natural to conclude, and we will see shortly that this is indeed the case, that the corresponding energy-momentum tensor is also rapidityindependent.
The most general energy-momentum tensor for a medium produced in a highenergy collision of two large nuclei with the rapidity-independent matter distribution can be parametrized in terms of the energy density and transverse and longitudinal pressures P T , P L as (2.8) (We employ translational invariance in the transverse plane due to the nuclei being very large.) At mid-rapidity (z = 0) the tensor (2.8) looks like in the t, x, y, z coordinates. It is also useful to note that the energy-momentum conservation The energy-momentum tensor in the massless ϕ 4 theory with the Lagrangian density is given by With the two-point correlation function decaying at late times as shown in Eq. (2.7), it is natural to conclude that the four-point correlation function decays twice as fast and is, therefore, negligible in Eq. (2.13) taken at late times. We therefore arrive at (2.14) Here the subscripts 1 and 2 denote derivatives with respect to x 1 and x 2 respectively. Substituting the Green function from Eq. (2.5) we find (cf. [53,57], the J 0 ↔ J 1 difference in p T is inconsequential at late times) for the classical medium at τ Q s 1. Applying the late-time asymptotics to Eqs. (2.15) we arrive at This is a free-streaming anisotropic medium, with zero longitudinal pressure and a non-zero transverse pressure. In the full MV model, the classical gluon fields produced in a nuclear collision also lead to the free-streaming asymptotics of Eq. (2.16). Note also that due to Eq. (2.11) the ∼ 1/τ scaling corresponds to P L = 0 and, hence, to free streaming. A more isotropic medium with non-zero P L > 0 would have the energy density falling off faster than 1/τ . The isotropization problem in heavy ion collisions can be formulated as follows: can one find (presumably quantum) corrections to the classical field correlator (2.3) modifying the free-streaming result (2.16) at late times in such a way as to generate a non-zero p L > 0, or, equivalently, give us the net energy density that decreases faster than 1/τ ?

Boltzmann equation prediction for the two-point function after a single rescattering
Kinetic theory appears to give a positive answer to the above isotropization question.
In [55] we show that using the classical particle distribution f cl = f (0) resulting from the correlator (2.1), inserting it in the collision term of the Boltzmann equation, and solving the resulting equation for the new distribution f (1) gives the following components of the energy-momentum tensor: Putting the strong coupling to zero in Eqs. (2.17), α s = 0, we recover the leadingorder classical free-streaming result (2.16) (with some coefficient A (0) ). The corrections due to the single iteration of the Boltzmann collision term enter at the order α 2 s in Eqs. (2.17) with some coefficients A (1) and B (1) . (Note that τ 0 is the lower cutoff for the applicability times of kinetic theory: it is applicable for τ > τ 0 .) Since the pressure components in kinetic theory must be positive, we see that B (1) > 0. Therefore, the B (1) -corrections in (2.17) appear to lead to P L > 0 and simultaneously make the energy density decrease faster than 1/τ . We see that kinetic theory predicts that a single 2 → 2 rescattering of two classical fields would generate the first isotropization correction we are after. The remaining question is whether such correction arises in the full field-theoretical calculation.
The perturbative solution of Boltzmann equation from [55] is quite general and is valid for the 2 → 2 collision term calculated in any theory. While the discussion in [55] concentrated on QCD with the expansion in Eqs. (2.17) being in powers of α s , it can be easily modified to apply to the scalar ϕ 4 theory in question here by simply replacing α s → λ everywhere. (All the appropriate constants and scattering amplitudes would be modified as well, but this would not matter so much for us since we are interested mainly in the form of the τ -dependence of the energy-momentum tensor.) Equations (2.17) can be rewritten as (1) and possibly even τ 0 are different in Eqs. (2.17) and (2.18).) Once again we observe that the kinetic theory has a specific prediction for the outcome of the single 2 → 2 rescattering of two classical fields. 1 In the remainder of this paper we will verify the results (2.18) by an explicit diagrammatic calculation.

Two-point correlation function with a single rescattering: Full diagrammatic calculation in momentum space
In this Section our aim is to calculate the late-times τ 1 , τ 2 asymptotics of the correlation function G 22 (x 1 , x 2 ) due to a single 2 → 2 rescattering involving two classical fields. The diagrams we want to calculate are shown in Fig. 3. The diagrams are labeled I, II, II', III, and III' with II' and III' obtainable from II and III by replacing x 1 ↔ x 2 in them. In the kinetic theory language, diagrams I, II, and II' correspond to the gain term in the (collision term of the) Boltzmann equation, while diagrams III and III' correspond to the loss term. Below we will first calculate diagrams I, II, and II' together, and then calculate diagrams III and III'.

Diagrams I, II, and II'
Let us start with diagram I. In momentum space one has with the leading-order classical correlation function given by (2.3). Substituting (2.3) into (3.1) and integrating out all the delta-functions (except for one) yields with k 2 = k − k 1 − k 3 and the retarded scalar Green function In coordinate space one has From now on it is implicitly understood that In arriving at Eq. (3.4) we have defined Equation (3.4) can be written more compactly by defining We obtain While the exact evaluation of I 3 appears to be rather involved, we can obtain its late-time asymptotics. To do this we start with the full expression for I 3 , and integrate over k − to obtain , in a form reminiscent of the light-cone perturbation theory (LCPT) [59]. Now let us apply the large-τ limit to Eq. (3.10). Note that since the limit of a product is equal to the product of the limits, we can first take the large-x + limit of the last fraction in Eq. (3.10). To do so we observe that in the distribution sense. Indeed this is true for any function h(a) of real variable a decomposable into a Fourier integral, To see this we simply point out that which is identical to the same convolution of e i a ξ with the right-hand side of Eq. (3.11), (3.14) Applying Eq. (3.11) to the last fraction in Eq. (3.10) after neglecting all the i 's in the exponent we obtain where by k − we now imply its on-shell value, Equation (3.15) easily simplifies into which, with the help of Eq. (3.6) becomes Here again k − is given by Eq. (3.16).
Note that here and throughout the paper, when writing τ → ∞, we mean large but finite proper time τ , that is τ p T 1 with p T being any of the transverse momenta in the problem.
For now we leave the diagram I as evaluated in Eq. (3.19) and turn our attention to diagrams II and II'. Employing Eq. (2.3) and integrating over all delta-functions except for one in diagrams II and II' gives in momentum space In coordinate space we have with k − given by Eq. (3.16). To study the late-time asymptotics of diagrams II and II' we employ Eq. (3.18). This gives Adding diagrams I, II, and II' given by Eqs. (3.19) and Eq. (3.22) together we arrive at I+II+II' The late-time asymptotics of diagrams I, II and II' given by Eq. (3.23) is dominated by the saddle points of the k + and k + integrals, unless the integrands have other singularities which may prevent deforming the integration contours into the steepest descent contours. Regardless of the analytic structure of the integrands, we can argue that late-time asymptotics at x + 1 = x + 2 , x − 1 = x − 2 (as is needed for calculation of expectation values of local operators, e.g. of the energy-momentum tensor) is dominated by the regions of integration where k + = −k + : indeed, the dominant values of k + and k + have to be such that the two oscillating exponentials in Eq. (3.23), , giving no oscillations in the end. (In fact, functions that oscillate rapidly at late times can be simply neglected in determining the asymptotics.) Therefore, expecting k + = −k + , we can put Sign(k + ) = −Sign(k + ) (3.25) in Eq. (3.23). This leads to and the sum of the diagrams I, II, and II' at late times becomes I+II+II' Such partial cancellation between the diagrams I, II, and II' was seen before in the framework of kinetic theory [25]. Further evaluation of the expression (3.27) appears to be impossible without an explicit expression for I 1 . The corresponding calculation is carried out in Appendix A resulting in where the branch cut of the square root is chosen along the negative imaginary axis for the later convenience. Note also that the sum k T + p T and the difference k T − p T involve the magnitudes of the transverse momenta, and are not a sum or a difference of vectors.
Once again let us point out that to obtain the late-time asymptotics of Eq. (3.27) we need to try to deform the contours of the k + and k + integrals into the steepest descent shape, in order to perform the saddle point approximation. This contour deformation may be affected by the presence of singularities in the k + and k + complex planes. For definitiveness, consider the k + integral. For all other momenta in Eq. (3.27) fixed, it has an essential singularity at k + = −i and branch cuts due to The steepest descent contour for the k + integral is shown in the left panel of Fig. 4. The saddle points of the k + integral are given by k + = ±k + sp with (3.29) They correspond to points (±1, 0) in both panels of Fig. 4. In the case of no singularities in the complex k + , it is clear that one can easily deform the k + integration from running along the real axis to the steepest descent curve in the left panel of Fig. 4. The right panel of Fig. 4 illustrates what happens if one tries to deform the real axis contour into the steepest descent one in the presence of singularities in the complex k + plane. In that plot we explicitly show the essential singularity at k + = −i in the integrand of Eq. (3.27): one can see that it does not interfere with the integration contour deformation into the steepest descent shape because it lies outside the contour, and, when the contour approaches this singularity, it does so along the positive imaginary axis near the origin, k + ≈ +i . This region of integration is exponentially suppressed as one can see by plugging k + ≈ +i into the first exponential of Eq. (3.27). Hence we do not need to worry about the singularity at k + = −i .
In contrast, the branch cuts may possibly interfere with the contour deformation: this is illustrated in the right panel of Fig. 4 by a sample vertical branch cut with the branch point on the real axis. This branch cut is not an accurate representation of the branch cuts of the k + integrand, and is shown here as a toy model to illustrate the possibility that in deforming the integration contour one may have to wrap the contour around a part of the branch cut as well. The corresponding sample integral would look like (cf. Eqs. (3.27) and (3.28)) with k + br − i the branch point near the real axis, as shown in the right panel of Fig. 4. (Note again the the branch cut of the square root is chosen along the negative imaginary axis.) Writing with some real variable y we can approximate the contribution to the integral in Eq. (3.30) from the part of the contour wrapped around the branch cut by (3.32) In the process we have assumed that the branch cut section of the contour is long enough for the upper limit of the y integral in Eq. (3.32) to be replaced by infinity: this is only valid if k + br is sufficiently far from the saddle point k + sp (and for |k + br | > |k + sp |). More specifically, we need |k + br − k + sp | 1/x − 1 ∼ 1/τ 1 . Since 1/τ 1 is small for large times τ 1 , this assumption is justified in most cases.
What we learn from the sample integration in Eqs. (3.30) and (3.32) is that the branch cut contribution is dominated by the branch point, k + = k + br , with a small region y ∈ [0, ∼ 1/τ 1 ] near the branch point contributing dominantly to (this part of) the integral.
We are now ready to tackle the full k + and k + integrals in Eq. (3.27). First we need to identify the branch points and branch cuts of the k + and k + integrals. Starting with the k + integral we see that its branch cuts originate in I 1 (k + − k + 3 , k − − k − 3 , k 1T , k 2T ) as follows from Eq. (3.28). Defining we conclude that there are four branch points given by The sign in (k 1T ± k 2T ) is either a plus or a minus simultaneously inside the square root and outside: one cannot have k 1T +k 2T in one place in the expression and k 1T −k 2T in another.) These branch points can be real or complex. In the latter case of complex ξ the branch point may contribute to the integral only if |Re ξ| > |k + sp |, as follows from the contour in Fig. 4. In such case one can easily show that the contribution of the complex-valued branch point in the lower k + half-plane is exponentially suppressed. Since for positive x − 1 that we are interested in one needs to close the k + integration contour in the lower half-plane, we conclude that we can discard the contributions of the complex-valued branch points.
We are left with the case of real-values branch points ξ. Due to the presence of θ(−k + k + 3 ) in Eq. (3.27), only the negative real ξ can contribute (otherwise the k + integration never approaches the branch point for it to contribute). A quick analysis of Eq. (3.34) shows that only two solutions can be real and negative. Let us denote them ξ 1 and ξ 2 , such that with ξ 2 < ξ 1 < 0. The corresponding branch point in the k + plane are given by k + = ξ 1 k + 3 and k + = ξ 2 k + 3 . Note that the branch points of the k + integral are given by the branch points of I 1 (k + + k + 3 , k − + k − 3 , k 1T , k 2T ) resulting in k + = −ξ 1 k + 3 and k + = −ξ 2 k + 3 . Here we will consider the case of k + 3 < 0 (with k + > 0, k + < 0): the k + 3 > 0 case can be done by analogy. The branch cuts of the k + and k + integral in the k + 3 < 0 case are shown in Fig. 5.
Remembering from the example above that the integrals around our branch cuts are dominated by the branch points, and invoking the argument that the late-time asymptotics is dominated by k + = −k + to avoid rapid oscillations (which would make the function practically zero), we conclude that if the k + integral is wrapped around the branch cut originating at, say, the ξ 1 k + 3 branch point, the k + integral must be wrapped around a branch cut originating at the −ξ 1 k + 3 branch point. However, as one can see from Fig. 5, while the branch cut starting at ξ 1 k + 3 in the k + complex plane lies in the lower half-plane, the branch cut starting at −ξ 1 k + 3 in the k + complex plane lies in the upper half-plane. For x − 1 , x − 2 > 0 that we are interested in one needs to close the k + and k + contours in the lower half-plane: hence one cannot simultaneously pick up the contributions of the branch cuts originating at ξ 1 k + 3 and −ξ 1 k + 3 (or at ξ 2 k + 3 and −ξ 2 k + 3 ) in the k + and k + integrals. 2 We thus conclude that both the k + and k + integrals can not be dominated by branch cuts.
For a given value of k + 3 , and for x ± 1 ≈ x ± 2 , the branch points are either near |k + br − k + sp | < ∼ 1/τ or far |k + br − k + sp | 1/τ from the saddle point. The 'near' region is small at large τ and its contribution is suppressed by an extra power of 1/τ . Hence the leading contribution to the integral in Eq. (3.27) comes from the 'far' region of |k + br − k + sp | 1/τ . In this region branch cuts are clearly separated from the saddle points: if we do not want to have rapidly oscillating exponentials, we can either have the contribution of saddle points in both the k + and k + integrals or of the branch points in both the k + and k + integrals, but due to the separation of branch cuts and saddle points we cannot have a mix where the branch point contributes in one and the saddle point contributes in the other. Since we have just eliminated the contribution of the branch points in both the k + and k + integrals, we are left with the contributions of saddle points only. The saddle point k + -integral gives We conclude that I+II+II' consistent with free streaming. 3 Note that the k + 3 integral in Eq. (3.38) is finite, as can be checked explicitly.

Diagrams III and III'
Our aim now is to calculate diagrams III and III' from Fig. 3. Following the steps we made for other diagrams, we write Substituting the leading-order correlators from Eq. (2.3), integrating out all but one delta-function, and also integrating over k + and k − explicitly yields Integrating over k − and employing the limit from Eq. (3.11) one arrives at where we have defined q ± ≡ k ± 1 + k ± 3 along with Consider the integral over q + and q − in Eq. (3.42): This object is boost-invariant. It is a function of several transverse momenta, and of only one longitudinal four-vector component -of k + . A boost-invariant object cannot depend on only one k + : hence, it must be independent of k + , 4 that is, The exact form of J(k T , k 1T , k 2T , k 3T ) is not important for the late-time asymptotics since now we can integrate over k + exactly, obtaining We conclude that Diagrams III+III' again consistent with free streaming.

Energy-Momentum Tensor
To cross-check our results let us calculate the longitudinal pressure at z = 0 (in the massless ϕ 4 theory at hand): At the leading saddle point k + sp = ± k T √ 2 e η 1 , k + sp = ∓ k T √ 2 e η 2 and at mid-rapidity (z = 0) the square brackets become (We have also put k = −k.) We see that the leading saddle point contribution gives the longitudinal pressure at late time, as characteristic of free streaming. For completeness, let us calculate the energy density: The saddle points now give (at z = 0) such that the energy density is again in agreement with the free-streaming behavior (2.16). Finally, the transverse pressure is Again, at the saddle points at mid-rapidity we get 56) and the transverse pressure is in complete agreement with the free-streaming expectations (2.16).

Two-point correlation function with a single rescattering: Full diagrammatic calculation in the Wigner representation
In this Section we calculate the late-time asymptotics of the correlation function G 22 (X, p) in the Wigner representation. We start with the calculation of the classical two-point correlation function given by the diagram in Fig. 2. Then, we calculate G 22 (X, p) due to a single 2 → 2 rescattering using the diagrams in Fig. 3.

The asymptotic expansion of the classical correlation function in
with k T = 0. One can integrate out k − by closing the integration contour downward and picking up the residues of the poles in the lower half-plane. There are two poles in the lower half-plane and their residues yield with the understanding that the k + integration contour is located infinitesimally above the real axis.
In this paper we are only interested in the asymptotic behavior of G 22 (X, p) at large τ . For this goal, we evaluate the two terms on the right-hand side of (4.3) separately. The integrand of each term has a steepest descent path passing through two saddle points respectively located at and k + = −2p + ± p T 2X + X − for the second term. (4.5) As illustrated in Fig. 6, we deform the integration contours for these two integrals so that they continue along the steepest descent paths. Each term has two other poles at Here, one only needs to consider the case with p 2 2p + p − > 0 since the residues of these poles will be exponentially suppressed at large τ otherwise. The contour deformation may involve passing through these poles. If this is the case, one needs to pick up  Figure 6. Deformation of integration contours for the 1 st (top) and 2 nd (bottom) terms on the right-hand side of (4.3). In each figure, the original integration contour is a horizontal line infinitesimally above the real axis. It is deformed so that it continues along the steepest descent path. The contributions from integrating over l 1 and l 2 cancel with each other and the contributions from the arcs, C 1 , C 2 and C 3 , vanish as their radius goes to ∞. the contribution from their residues. By integrating over the deformed contours one can easily pick up the contributions from the residues of the poles (4.6) and from the steepest descent paths.
The residues of the poles yield a contribution with an amplitude of O(τ 0 ) at large τ . It is easy to see that the integrand on the right hand side of (4.3) as a whole is not singular at the points in (4.6). That is, the residues of the two terms at the same pole cancel with each other. Therefore, we only need to take into account the case when the deformed contours for the two terms pass through different poles: k + = 2p + p 2 2p + p − for the first term and k + = −2p + p 2 2p + p − for the second term. As illustrated in Fig. 6, for p + > 0 this requires which is equivalent to Similarly, for p + < 0 this requires which is also equivalent to (4.8). By combining the results for both p + > 0 and p + < 0, we obtain the following contribution to G LO 22 (X, p) from the residues of the poles Poles = cos 2 As shown in Fig. 7, the above result agrees very well with our numerical result of G LO 22 (X, p) at large τ . The integration over the steepest descent paths gives a contribution with an amplitude of O(τ − 1 2 ) at large τ . Such a contribution comes from the integration over a small region in the vicinity of each saddle point in (4.4) and (4.5). One only needs to expand the exponent of the exponential function in a Taylor series around each saddle point up to the second order and replace k + in the rest part of the integrand by the saddle point value. By doing this and integrating out k + , we get As shown in the figure on the top of Fig. 7, the above result is indeed needed for us to understand the τ dependence of G LO 22 (X, p) at an intermediate large time except when the saddle points coincide with poles. This coincidence explains the divergence between the exact numerical results and the dash-dotted line in the top panel of Fig. 7. Our exact numerical results in Fig. 7 show that G LO 22 (X, p) is well behaved in these regions. Hence, these regions should not be important for calculating any observables involving G 22 (X, p) and we do not need to construct a separate analytical expression for these regions. As one can see from the lower panel of Fig. 7, these regions become progressively less important at later times.
In Ref. [55] we find that the classical gluon two-point function is we obtain from Eq. (4.10) (4.15) for our scalar correlator. In order to get the correct coefficient of the above equation, we write (4.16) and fix c θ by integrating out p ± and matching it to that by integrating (4.10) over p ± . It is convenient to define ν ≡ p 2 = 2p + p − − p 2 T and y = 1 2 ln p + p − . In terms of these two variables, Eq. (4.10) reduces to, for p + p − > 0, (4.17) and, for p + p − < 0, with only Eq. (4.17) contributing to n 1 and only Eq. (4.18) contributing to n 2 . After changing variables toŝ = sinh(y − η) andν = p 2 /p T , we have We only need to obtain the large τ asymptotics. At large τ , the predominant contribution to the above integral comes from the regionν ≈ 0 andŝ ≈ 0. So one can not take either τν or τŝ as the large expansion parameter. We do the following trick: we rotate the integration contour of, sayŝ, to go along the positive imaginary axis. Then the [0, i] region of integration does not contribute to the real part of the integral. We get Similarly, by changing variables toŝ = cosh(y − η) andν = −p 2 /p T , we have One can simply discard n 2 since it is a highly oscillatory function at large τ . From n 1 we obtain c θ = 1/2 and where with k = 0 and k 2 = p−k 1 −k 3 . The integration contours for k + and k − are understood to be located infinitesimally above the real axis of the k ± plane and we shall drop all the i 's in both G R and I 1 . In order to calculate the late-time asymptotics, we first need to identify all the critical points including poles, branch points and saddle points. For k − , the integrand on the right hand of (4.24) has two poles given by (4.2) and 4 branch points at and the solid lines in Fig. 8. At late times, the poles give a contribution proportional to τ 0 while the branch cuts yield a contribution proportional to τ − 1 2 . Hence, we only keep the contribution from the residues of the poles. As a result, we get (4.28) For k + , the integrand on the right hand side of (4.28) has poles given by (4.6). It also has saddle points and branch cuts. Since they only give a contribution proportional to τ − 1 2 , we only need to keep the contribution from the residues of the poles. The calculation is straightforward and similar to the leading-order case of Sec. 4.1. We obtain where we have used the identity to shorten the expression by taking the real part. Based on the same calculation as for (4.23) we obtain  Figure 9. Relative positions of poles with respect to the steepest descent path. Here, we assume that p + > 0 and that the two poles are located at k + = ±2p + p 2 2p + p − . These 3 figures show all the possible cases when one needs to keep the residues of the poles in the k + integration in (4.33).

From (3.20) we have in Wigner representation
with k = 0 and k 2 = p − k 1 − k 3 . By integrating out k − we obtain with k − 1p and k − 2p given in (4.2). By using Eq. (4.30) one can show that the second term on the right hand side of (4.33) gives a contribution which is the complex conjugate of the first term.
At the end we have (4.41) That is, at late times Combining (4.37) and (4.42) we conclude once again that the late-time asymptotics of the rescattering diagrams in Fig. 3, calculated this time in Wigner representation, is consistent with free streaming.

Conclusions and Outlook
In the first paper [55] of this 'duplex' we have outlined the way to apply the Schwinger-Keldysh formalism to ultrarelativistic heavy ion collisions, thus setting the stage for the calculation of time-dependent observables, such as the energy-momentum tensor, in the collisions. As the first application of this technique, in [55] we have tried to re-derive the Boltzmann equation for the medium produced in heavy ion collisions by considering a single 2 → 2 rescattering correction to the classical gluon fields of the MV model. We have employed the "on-shell" approximation for the propagators that one usually employs in deriving the kinetic theory. The result was dependent on the time-ordering assumptions outlined above in Sec. 1: kinetic theory emerged under the assumption (i), while free-streaming was obtained if assumption (ii) was employed.
In this paper we have used the formalism set up in [55] to re-do the calculation of the 2 → 2 rescattering correction to the classical fields without using the "on-shell" approximation for the propagators and without assuming a specific time-ordering of the interaction time versus the measurement time. Performing the calculation twice, both in momentum space (Sec. 3) and in the Wigner representation (Sec. 4), we have arrived at the results consistent with free streaming (2.16). We have thus found no evidence for the applicability of the kinetic theory (employing the Boltzmann equation taken with the full collision term) to the perturbative description of heavy ion collisions.
Our perturbative calculations of the energy-momentum tensor are consistent with the conclusion reached in [53], where it was argued that the energy density of the produced weakly-coupled medium is given by (τ, η, b) at any order in perturbation theory. In the scenario advocated in [53], higher-order perturbative corrections would only modify the gluon multiplicity distribution dN/d 2 k T dη d 2 b T and the saturation scale Q s , leaving the τ -dependence of Eq. (5.1) unchanged.
In the future we hope that the formalism we have presented in [55] will be useful for calculations of time-dependent heavy ion observables. In particular it can be used to cross-check (and possibly challenge) the conclusion in (5.1) by explicit perturbative calculations. Our own cross-check presented above did not show any deviations from (5.1) and disagreed with kinetic theory. Perhaps other thermalization scenarios may fare better in challenging Eq. (5.1). Indeed calculations of higher-order perturbative corrections to the classical gluon field contributions to heavy ion observables appear to be very complicated. In our minds, however, such calculations would present a necessary theoretical test for any thermalization scenario. If, for instance, a given thermalization proposal claims to resum a certain τ -dependent parameter to all orders, then this parameter should manifest itself in some lower-order perturbative calculation. In other words, one needs to prove that the resummation parameter exists. If the explicit calculation fails to confirm that the resummation parameter exists, the thermalization scenario in question should be discarded. If the equation (5.1) does turn out to be exact in the perturbation theory, then calculation of higher-order corrections would improve our knowledge of gluon multiplicity and energy density, which would also be very useful. We hope the formalism and calculations presented in [55] and in this paper lay the groundwork for future cross-checks of thermalization scenarios and will facilitate calculations of heavy-ion observables in the years to come.

A I 1 Calculation
The goal of this Appendix is to calculate Taking G R from Eq. (3.3) and integrating (A.1) over k − we get 3) The poles of the integrand are given by Using k + 1,2 we rewrite Eq. (A.3) as , (A. 6) where the coefficients r 1 , r 2 have to be determined by matching the linear in terms in the denominators of (A.6) and (A.3) at the poles k + = k + 1 and k + = k + 2 respectively. This gives Multiplying Eqs. (A.7) after a little algebra we arrive at (A.8) We are now ready to integrate Eq. (A.6) over k + . To do so we need to consider two cases: • Case I: D < 0. In this case k + 1 and k + 2 have non-zero (and non-infinitesimal) imaginary parts, and the i terms along with the values of r 1 and r 2 in Eq. (A.6) are not important. Integrating Eq. (A.6) over k + we obtain We conclude that (A.12) By combining the above cases I and II, we have In order to perform integrations using the complex plane, it is desirable to rewrite I 1 without θ-functions, since the latter are hard to analytically continue into the whole complex plane. This is achieved by re-writing Eq. (A.13) as 2(q + +i ) 1 2 . (A.14) The i regulators are needed for the "standard" branch cut of the square root running along the negative real axis. If we choose the branch cut of the square root to run along the negative imaginary axis, one can neglect these i 's, arriving at the Eq. (3.28) in the main text.