Kinetic theory of a longitudinally expanding system of scalar particles

A simple kinematical argument suggests that the classical approximation may be inadequate to describe the evolution of a system with an anisotropic particle distribution. In order to verify this quantitatively, we study the Boltzmann equation for a longitudinally expanding system of scalar particles interacting with a $\phi^4$ coupling, that mimics the kinematics of a heavy ion collision at very high energy. We consider only elastic $2\to 2$ scatterings, and we allow the formation of a Bose-Einstein condensate in overpopulated situations by solving the coupled equations for the particle distribution and the particle density in the zero mode. For generic CGC-like initial conditions with a large occupation number and a moderate coupling, the solutions of the full Boltzmann equation do not follow a classical attractor behavior.

1 Introduction and motivation 1

.1 Context
A long standing problem in the theoretical study of heavy ion collisions is the time evolution of the pressure tensor and its isotropization [1,2,3,4,5,6,7,8,9,10,11,12]. This question is closely related to the use of hydrodynamics in the modeling of heavy ion collisions. Although the complete isotropy of the pressure tensor is not necessary [13] for the validity of hydrodynamical descriptions, the 1 arXiv:1506.05580v2 [hep-ph] 29 Jun 2015 ratio of longitudinal to transverse pressure, P L /P T , should increase with time for a smooth matching between the pre-hydro model and hydrodynamics.
In most models with boost invariant initial conditions, the longitudinal pressure is negative immediately after the collision [14,15] (in the Color Glass Condensate framework [16,17,18,19,20], this can be understood as a consequence of large longitudinal chromo-electric and chromo-magnetic fields). On timescales of the order of a few Q −1 s (Q s is the saturation momentum), the longitudinal pressure rises, and reaches a positive value of comparable magnitude to the transverse pressure. However, it is unclear whether the scatterings are strong enough to sustain this mild anisotropy, or whether the longitudinal expansion wins and causes the anisotropy to increase.
This regime may be addressed by various tools. Indeed, it corresponds to a period where the gluon occupation number is still large compared to 1, but small compared to the inverse coupling g 2 so that a quasi-particle picture may be valid. This means that the system could in principle be described either in terms of fields or in terms of particles. Since the occupation number is large, it is tempting to treat the system as purely classical. In a description in terms of fields, this amounts to considering classical solutions of the field equations of motion, averaged over a Gaussian ensemble of initial conditions whose variance is proportional to the initial occupation number. In a description in terms of particles, i.e. kinetic theory, this amounts to keeping in the collision integral of the Boltzmann equation only the terms that have the highest degree in the occupation number [21,22,23].

Classical attractor scenario
The field theory version of this classical approximation has been implemented recently [24] for scalar fields with longitudinal expansion, and it leads to a decrease of the ratio P L /P T , like the power τ −2/3 of the proper time. This behavior seems universal; it has been observed for a wide range of initial conditions, and both for Yang-Mills theory and scalar field theories with a point-like interaction (such as a φ 4 interaction) [1,2,3]. Based on this observation, it was conjectured that the time evolution of P L /P T can be decomposed in three stages, nicely summarized in Figure 3 of ref. [25] : i. A transient stage that depends on the details of the initial condition; ii. A universal scaling stage, called "classical attractor," during which the dynamics is purely classical and the ratio P L /P T decreases as τ −2/3 ; iii. A final stage where the occupation number has become of order 1, and where quantum corrections are important. The isotropization of the pressure tensor may happen during this final stage.

Classical approximation in anisotropic systems
However, a possible difficulty with the classical approximation is that the occupation number cannot be large uniformly at all momenta, which may make this approximation unreliable since the dynamics integrates over all the momentum modes. In particular, this may be the case in anisotropic systems, as we shall explain now in a kinetic theory framework. Suppose for discussion that the particle distribution has become extremely anisotropic, with a narrow support in p z , This situation occurs in the pre-equilibrium stage of heavy ion collisions, that can be described at leading order by a rapidity independent classical color field. In a non expanding system 1 , we expect that 2 → 2 scatterings will kick particles out of the transverse plane and restore the isotropy of the distribution. The Boltzmann equation that describes the evolution of f (p ⊥ , p z ) should be able to capture this isotropization. If we track the momentum p 4 (see Figure 1), the Boltzmann equation can be sketched as follows: We have separated the purely classical terms (cubic in f ) from the subleading terms that are quadratic in f . The classical approximation amounts to dropping the quadratic terms. The usual justification of this approximation is to say that the cubic terms are much larger than the quadratic ones when f 1. However, as we shall see now, this counting is too naive when the support of f is anisotropic.
If p 1 and p 2 are in the transverse plane, then p z3 + p z4 = 0. Therefore, if we request that the particle 4 is produced outside of the transverse plane, then we have f 3 = f 4 = 0. Because of this, many pieces of the collision term (underlined in eq. (2)) are zero. In particular, all the classical terms vanish. The physical interpretation is clear: the f 3 terms correspond to stimulated emission, which cannot happen when both final particles lie in an empty region of phase space.
The only non-zero term is the one in f 1 f 2 , but one must go beyond the classical approximation in order to capture it. This is true no matter how large the distribution f (p) is inside its support, i.e. even if the classicality condition f (p) 1 is satisfied there. Moreover, this argument can be trivially generalized to any n → n scattering process in kinetic theory 2 . Therefore, when the particle distribution is anisotropic, the classical approximation artificially suppresses out-of-plane scatterings at large angle, possibly resulting in wrong conclusions regarding isotropization. Now let us relax the δ function assumption, and consider the case where the range of angles with large occupancy is finite but narrow. In particular, suppose that for momenta p ∼ Q, the particles mostly reside with |p z | < δQ, where δ 1 describes how anisotropic the momentum distribution is. We are interested in scattering processes which move particles out of this highly-occupied region. So consider a 2 ↔ 2 scattering process, where the final momentum p 4 lies outside this highly-occupied region, so f 4 is small. Within the classical approximation, Eq. (2) is dominated by the f 1 f 2 f 3 term. But for all three of these occupancies to be large, we need |p 1z |, |p 2z |, |p 3z | < δQ. Since p z1 + p z2 = p z3 + p z4 , this ensures that |p z4 | < 3δQ, that is, the final state particle produced in a classical scattering must still have quite a small |p z | value. We will refer to these as classical scatterings. A scattering with |p 3z | ∼ Q and |p z4 | ∼ Q will always be suppressed by a small occupancy in the classical approximation, so these scatterings can be neglected classically, and only occur because of the quantum f 1 f 2 term in Eq. (2). We will call these quantum scatterings. Now let us see whether scatterings with |p z3 |, p z4 | ∼ Q may still be important, even if the occupancies f (|p z | ∼ δQ) 1 are large. First of all, while the integrand in Eq. (2) is small for quantum scatterings, the integration phase space is much larger for these processes. Specifically, since |p z1 |, |p z2 | ∼ δQ in all cases, and because of the momentum-conserving delta function, the phase space for "quantum" scatterings is larger than that for "classical" scatterings by a factor of ∼ 1/δ. Therefore, even if the occupancy in the highly-occupied region is f ∼ 1/δ, order-1 of the scatterings are "quantum." That is, order-1 of the scatterings involve quantum effects when the angle-averaged occupancy 2 The classical approximation of a n → n collision term contains only terms of degree f n+n −1 . Because of longitudinal momentum conservation, if one of the momenta points out of the transverse plane, then there must be at least one out-of-plane momentum. Therefore, of the n + n − 1 distribution functions, at least one is zero. below p ∼ Q is order 1. This argument may not apply in gauge theories, where the matrix element is strongly enhanced for small-angle processes, but it should be valid in a scalar theory where the matrix element is isotropic.
In addition, for many purposes, an individual scattering resulting in |p z4 | ∼ Q is much more important than a scattering resulting in |p z4 | ∼ 3δQ. A particle's contribution to the longitudinal pressure involves p 2 z /p 2 . The particle produced in a "quantum" scattering, with |p z | ∼ Q, contributes more to the longitudinal pressure than the classically-scattered |p z | ∼ δQ particle, by a factor of ∼ 1/δ 2 . Therefore, in terms of generating longitudinal pressure, each "quantum" scattering contributes a more important effect, by a factor of 1/δ 2 . Combining this factor with the larger phase space for quantum scattering, we find that, for the purposes of understanding the longitudinal pressure, the classical approximation already starts to fail when f (p z ∼ 0, |p| ∼ Q) ∼ δ −3 , which is when the angle-averaged occupancy is 1/δ 2 1. This suggests that the classical approximation may lead to missing some contributions that are important for isotropization in theories where the crosssection is dominated by large-angle scatterings, such as a φ 4 scalar theory. When the classical approximation is used in this theory, it has been observed in ref. [24] that the ratio P L /P T decreases like τ −2/3 at late times for a longitudinally expanding system (dotted curve in Figure 2). Given the above argument, the f 2 terms that are neglected in the classical approximation could alter this behavior sooner than one might naively expect, ending the growth of anisotropy and giving a ratio P L /P T ∼ τ 0 at late times. If the f 2 terms are truly negligible over some extended period of time, then the full Boltzmann equation should lead to the red curve in Figure 2, in which the system spends some time stuck  into a classical attractor before eventually leaving it in order to isotropize (this red curve corresponds to the 3-stage scenario of Section 1.2). In contrast, if the f 2 terms are important from the start, one may get the behavior illustrated by the blue curve, where the power law τ −2/3 does not play any particular role in the evolution of P L /P T and a classical attractor would not exist.

Contents
In the rest of this paper, in order to assess quantitatively this issue, we solve the Boltzmann equation for an expanding system of scalar bosons with a φ 4 interaction, both for the full collision kernel and its classical approximation. In Section 2, we discuss the Boltzmann equation for a longitudinally expanding system and prepare the stage for its numerical resolution. We describe our algorithm in Section 3, and the numerical results are exposed in Section 4. Section 5 is devoted to a summary and concluding remarks. Some more technical material and digressions are presented in appendices. In appendix A, we present results on the isotropization of the pressure tensor in a non-expanding system, that corroborate the fact that the f 2 terms play an essential role. In appendix B, we derive an analytical expression for the azimuthal integrals of the 2 → 2 collision term. Additional details about our algorithm can be found in appendices C and D. In the appendix E, we present an alternate algorithm for solving the Boltzmann equation, based on the direct simulation Monte-Carlo (DSMC) method.

Boltzmann equation for expanding systems 2.1 Notation
In the following, we consider partially anisotropic particle distributions that have a residual axial symmetry around the z axis, (p ⊥ ≡ p ⊥ ). For simplicity, we do not write explicitly the space and time dependence of f . The function f depends smoothly on momentum, except possibly at p z = p ⊥ = 0 if there is a Bose-Einstein condensate (see Section 2.4).
We also assume that f is even in the longitudinal momentum p z The Lorentz covariant form of the Boltzmann equation reads where C p [f ] is the collision term (see the subsection 2.3). This form of the equation can then be specialized to any system of coordinates. We will focus on the collision term which arises at lowest order in the coupling, which for φ 4 theory is elastic 2 ↔ 2 scattering. For a system that expands in a boost invariant way in the longitudinal direction, the most appropriate system of coordinates is (τ, η, x ⊥ ) for the spacetime coordinate and (y, p ⊥ ) for the momentum of the particle, which are related to the usual Cartesian coordinates and momenta by The assumed boost invariance of the problem implies that the distribution f does not depend separately on η and y, but only on the difference y − η. Therefore, it is sufficient to derive the equation for the distribution at mid-rapidity, η = 0. For simplicity, we will further assume that the distribution is independent of the transverse position.

Free streaming term
In the system of coordinates defined in eqs. (6), the left-hand side of the Boltzmann equation can be rewritten as where we have defined In a boost invariant system, f depends on y and η only via the difference y − η and it is sufficient to consider only the distribution at mid-rapidity, η = 0, for which we have This formula implicitly assumes that the distribution f is given in terms of the transverse momentum p ⊥ and the rapidity y. Energy and momentum conservation will take a simpler form if we express it in terms of the energy ω p and longitudinal momentum p z , instead of p ⊥ and y. Since p z = M p sinh(y) and ω p = M p cosh(y), one gets 3 3 The more familiar form is obtained when one expresses f in terms of pz and p ⊥ .

Collision term
If we consider only 2 → 2 scatterings, the right-hand side of the Boltzmann equation contains integrals over the on-shell phase spaces of three particles. This 9-dimensional integral can be reduced to a 5-dimensional integral by using energy and momentum conservation. A further simplification results from our assumption that the distribution is invariant under rotations around the p z axis. The Boltzmann equation reads with In this equation, p denotes the integration over the invariant phase-space and F nc ({P i }) is the factor that contains the particle distribution 4 , (We use the abbreviation f i ≡ f (p i ).) The invariance under rotation around the p z axis can be used in order to perform analytically all the integrals over azimuthal angles (we use the same procedure as in the case of an O(3) rotational invariance, see refs. [26,27]). In order to achieve this, let us first write By combining this with the integrations over the transverse momenta, we get The integral over x ⊥ depends only on the transverse momenta {p ⊥i }, but not on the distribution function f . It can therefore be calculated once for all. In the appendix B, we obtain an explicit formula for this integral in terms of the Legendre elliptic K function. Defining we obtain In terms of this integral, the collision term reads +∞ m dω p2 dω p3 dω p4 I 4 ({p ⊥i }) dp z2 dp z3 dp z4 For the sake of brevity, we have not written explicitly the boundaries of the integration range on p z2 , p z3 , p z4 . They are given by Eq (20) reduces to a 4-dimensional integral after taking into account the two delta functions, which is doable numerically. Assuming that we encode the momenta with (ω p , p z ) there is no need to perform any interpolation when using the delta functions to eliminate two of the integration variables, which is useful for fulfilling with high accuracy the conservation of energy and particle number. After this reduction that eliminates the variables ω p2 and p z2 , eq. (20) reduces to Cp z dp z3 dp z4 At each time step, C p1 [f ] must be calculated for each value of p z1 and ω p1 . Therefore, if we discretize the variables ω and p z with N f and N z lattice points respectively, the computational cost for each time step scales as (N f N z ) 3 .
In the numerical implementation, the allowed energy and momentum ranges must be bounded. We will denote the maximum allowed energy as ω Λ and the maximum allowed longitudinal momentum as L, so that The integration domain C ω for ω p3 and ω p4 is the bounded domain shown in Figure 3. Likewise, if we assume that p z1 > 0 (since f is even in p z ), the integra- tion domain C pz for p z3 and p z4 is the domain shown in Figure 4. The peculiar shape of these domains comes from the fact that ω p2 and p z2 extracted from the delta functions must themselves have values which lie within the bounds (23). The variable p z can be positive or negative, and although f is even in p z , we must integrate over both p z3,4 ≥ 0 and p z3,4 ≤ 0 inside the collision term, because of the asymmetric shape of the integration domain. Note that there is also an extra constraint (not represented on these diagrams because it relates ω and p z ) that must be satisfied, for the transverse momentum to be defined.

Bose-Einstein condensation
If only elastic collisions are taken into account, a Bose-Einstein condensate (BEC) may appear in the system if the initial condition is overpopulated [4,28,29,30,31]. When this happens, the particle distribution has a singularity at zero momentum, that we can represent by 6 After this redefinition, f describes the smooth part of the occupation number and n c is the particle density in the condensate. The Boltzmann equation (11) needs to be revised to account for n c . Indeed, one can now have particlecondensate interactions. One can easily check that the 2 ↔ 2 processes cannot involve more than one particle from the condensate. The modified Boltzmann equation thus reads [f ] is the contribution from collisions between the particle of momentum p 1 that we are tracking and a particle from the condensate, while C 12↔c4

p1
[f ] stands for a collision between the particle p 1 that we are tracking and another particle of momentum p 2 to give a final state with a particle in the condensate and a particle of momentum p 4 , This term should be doubled, to account for the fact that the condensate particle can be the particle 3 or the particle 4. Following the same procedure 7 as in Section 2.3, we find 6 Our convention for the integral of a delta function over the positive real semi-axis is so that the integral of the second term with the measure d 3 p/(2π) 3 gives nc. 7 Here one can directly use the formula (73) for the integral of three Bessel functions.
where A(x, y, z) is the area of the triangle of edges x, y, z : Note that in these collision terms, the 2-dimensional integration domains of Figures 3 and 4 collapse to 1-dimensional domains (see the appendix C). The equation for the evolution of the particle density n c in the condensate can be obtained from eq. (22), by replacing the particle p 1 that we are tracking with the singular part of eq. (25), and then integrating over p 1 . By doing so we obtain the following equation for the condensate

Conservation laws
The Boltzmann equation fulfills several conservation laws, that play an important role in determining the form of the equilibrium particle distribution and also in assessing the accuracy of algorithms employed for solving the equation numerically. Each collision conserves energy and momentum. In addition, since we are only considering elastic scatterings in this paper, the collisions also conserve the number of particles.
The particle density is given by Note that this is a density of particles per unit of rapidity η. Since a given interval of η corresponds to a volume that expands linearly with the proper time τ , the conservation of the total number of particles implies that Likewise, the components of the energy-momentum tensor are given by (p µ ≡ (ω p , p ⊥ , p z ).) In a longitudinally expanding system, the conservation of energy and momentum, ∂ µ T µν = 0, becomes

13
where ≡ T 00 is the energy density and P L ≡ T 33 is the longitudinal pressure. The fact that the solutions of the Boltzmann equation satisfy the conservation equations (33) and (35) is a consequence of the delta function δ(ω p1 + ω p2 − ω p3 −ω p4 ) and of the symmetries of the collision term under various exchanges of the particles 1, 2, 3, 4. Namely, the integrand in the collision term is symmetric under the exchange of the initial state or final state particles : and antisymmetric if we swap the initial and final states: Therefore, any approximation scheme used in a numerical algorithm should aim at satisfying these properties with high accuracy. The easiest way to implement the delta functions without loss of accuracy is to use a lattice with a constant spacing in the variables ω p and p z . By doing this, one is guaranteed that the values of ω p and p z obtained by solving the constraints provided by the delta functions are also points on this grid. In addition, the quadrature formulas used for approximating the integrals in the collision term should lead to which are consequences of the symmetries (36) and (37). In other words, even if these symmetries are manifest in the collision kernel, one should be careful not to violate them with an improper choice of the quadrature weights. As we will see, our numerical scheme respects these symmetries, so that the conservation laws can be satisfied up to machine precision. It is also important to realize the role played in the conservation laws by the terms that appear on the left-hand side of the Boltzmann equation. For instance, these terms provide the term ( + P L )/τ in eq. (35), since we have However, this identity relies on a cancellation between the boundary terms that result from the integration by parts on ω p1 and p z1 . This must be kept in mind when discretizing the free streaming part of the Boltzmann equation, in order to avoid introducing violations of the conservation laws through these boundary terms. 14

Classical approximation
A central question in this paper is the interplay between the classical approximation and isotropization. Each computation will therefore be performed twice: i. With the full expression for the combination of distribution functions that appear in the equations (14), (29) and (31); ii. In the classical approximation where only the cubic terms in the particle distribution are kept. This entails the following changes: Eq. (14): Eq. (29): As will become clear in the description of our algorithm in the next section, we use fixed cutoffs on the energy ω p ≤ ω Λ and on the longitudinal momentum |p z | ≤ L. These cutoffs are not exactly the same as in the implementation of the classical approximation in classical lattice field theory, where one uses fixed cutoffs on the transverse momentum and on the Fourier conjugate ν of the rapidity. Indeed, ν ≈ p z τ , so that a fixed cutoff on p z roughly corresponds to a cutoff on ν that grows linearly with time. Fortunately, we do not expect physical effects from these cutoffs if they are taken high enough, since the occupancy typically falls off exponentially at large energy and momentum.
Note that there is a variant of the classical approximation defined in (41), in which each distribution function f is replaced by f + 1 2 . It is well known that this Ansatz provides the correct quadratic terms [21,22], accompanied by some spurious terms that are linear in the distribution function. This variant is known to suffer from a severe ultraviolet cutoff dependence, when the cutoff becomes large compared to the physical scales [32] (this property is closely related to the non-renormalizability of a variant of the classical approximation in quantum field theory, where one includes the zero point vacuum fluctuations [33]). For this reason, our algorithm for the Boltzmann equation in a longitudinally expanding system cannot be employed to study this alternate classical approximation, because of its fixed cutoff in p z , while the physical p z s decrease with time due to the expansion. In appendix A, where we consider the question of isotropization in a non-expanding system, we have also included this variant of the classical approximation (labeled "CSA" in Figures 15 and 16) to the comparison with the full calculation, and it appears that the quadratic terms in f included within this variant considerably improve the agreement with the full Boltzmann equation regarding isotropization. 15 3 Algorithm

Discretization
We adopt the following discretization for the longitudinal momentum and the energy: • The longitudinal momentum p z is taken in the range [−L, L], which we discretize into 2N z + 1 points (including the endpoints p z = ±L). The longitudinal step ∆p z and the discrete values are given by • The energy ω p is taken in the range [m + ∆ω, ω Λ ], which is discretized in N f points. The step ∆ω and the discrete values • The transverse momentum p ⊥ is then defined as If the inequality is not satisfied, then the pair (i, j) is excluded from the lattice.
• The particle distribution f (p) is encoded as a function of ω p and p z . In addition, the assumed parity of f in p z translates into The motivation for this choice is that, with uniformly spaced discrete energy values, a scattering from a pair of lattice points to a pair of lattice points will exactly represent energy conservation. Further, an integral over all momenta can be represented as a sum over lattice positions; no sampling or Monte-Carlo integration errors ever arise, only errors from the discretization procedure itself.

Collision term
Let us now present an algorithm that preserves all the symmetries of the collision kernel. This algorithm is a simple extension of the one used in ref. [32] (appendix B.2). For simplicity, we can assume that p z1 > 0, since the particle distribution is even in p z . The integrals over ω p and p z in the collision kernel are approximated by two 1-dimensional quadrature formulas. For the integral of a function F (ω p ), we use In our implementation, we have chosen the weights w f [i] as follows 8 Similarly, for integrations over the longitudinal momentum p z we use with the following weights If we denote (i 1 , j 1 ) the integers corresponding to the momentum p 1 , the expression (20) for the collision kernel in the absence of condensate can thus be approximated by In these sums, one should discard any term for which one of the pairs (i a , j a ) does not comply with the inequality (44). Eqs. (38) have been checked to hold with machine accuracy with this scheme. Similar formulas can be written for the terms that describe collisions involving a particle from the condensate, and they are given in appendix C.

Free streaming term
In the absence of collisions (e.g., in the limit g 2 → 0), the Boltzmann equation describes free streaming, a regime in which each particle moves on a straight line with a constant momentum. In the system of coordinates (τ, η), we are describing a slice in the rapidity variable η. This slice is progressively depleted of its particles with a non-zero p z , since they eventually escape, and the support of the particle distribution in p z therefore shrinks linearly with time. On our lattice representing discrete values of ω p and p z , the derivatives with respect to ω p and p z that appear on the left-hand side of the Boltzmann equation represent the fact that each particle systematically loses p z , and therefore energy ω p . The trajectory of a particle in (p z , ω p ) space, shown in Fig. 5, will be inward and downward. Discretizing the allowed p z , ω p values, this can be viewed as the particles hopping inward and downward. For instance, the particle number on the site (i, j > N z ) will move towards the sites 9 (i, j − 1) and (i − 1, j − 1). Meanwhile, particle number living on the site (i, j + 1) or (i + 1, j + 1) will flow onto the site (i, j). In contrast, a particle on the site (i, j < N z ) (i.e. p z < 0) can hop to (i, j +1) or (i−1, j +1), while a particle located at (i, j −1) or (i+1, j −1) can jump to (i, j). Once a particle reaches the line j = 0, it does not move from there. These moves are illustrated in Figure 5.
, it is sufficient to consider j ≥ 0. Let us denote the number density of particles at site (i, j) as The most general form for a discrete version of the collisionless Boltzmann equation can be written as 10 9 This choice is not unique. For instance, one could instead consider hops to (i, j − 1) or (i − 1, j), which would correspond to a different discretization of the derivatives ∂ω and ∂p z . We choose this discretization because the pz change is always larger than the ωp change, and because the point (i − 1, j − 1) almost always exists, while along the kinematic boundary the point (i − 1, j) generally does not. 10 The special case j = 0 is written explicitly in the eq. (92) of the appendix D.
where (α, β, γ) ij are coefficients that will be adjusted in order to satisfy all the conservation laws 11 . The non-condensed contribution to the particle density reads Similarly, we can write its contribution to the energy density and longitudinal pressure as follows In order to fully determine the unknown coefficients, we also need to consider the first moment of the distribution of longitudinal momenta, In the appendix D, we show that (α, β, γ) ij must have the following form 12 : Starting from (52) and replacing the local particle density by (51), we obtain the following discretization (for j > 0) for the free streaming equation 11 We can disregard the condensate in this subsection. Indeed, since the particles in the condensate have zero momentum, they play no role in free streaming, which simply causes the condensate number density to decay as τ −1 .
12 A limitation of our scheme is that it works only if the points (i = 1, j = ±1, ±2, · · · ) are not allowed by the condition ω 2 This can be fulfilled by choosing appropriately the lattice parameters (these points have been surrounded by a circle in Figure 5 -in the example of this figure, the above inequality is not satisfied). All the numerical results shown in this paper have been obtained with a lattice setup that satisfies this condition.
In the particular case j = 0, this equation reads One can also rewrite eq. (58) as For the internal points, where all the weights w f [i] and w z [j] are equal to one, this becomes which indeed reproduces the left-hand side of the Boltzmann equation (26) in the continuum limit.

CGC-like initial condition
We now solve the coupled equations (26) and (31) with the algorithm described in the previous section. We have used a moderately anisotropic initial condition that mimics the gluon distribution in the Color Glass Condensate at a proper time τ ∼ Q −1 s . It is characterized by a single momentum scale Q, below which most of the particles lie. The scale Q also sets the unit for all the other dimensionful quantities. To be more specific, our initial distribution at the initial time Qτ 0 = 1 is: with a large occupation below Q, f 0 = 100. According to the standard argument, such a large value of f 0 should ensure that the classicality condition is well satisfied, and that the cubic terms alone lead to good approximation of the full solution. We choose a coupling constant 13 g 4 = 50. The particle density in the condensate is initially n c,init = 10 −6 , and the mass m is taken to be m/Q = 0.1. Finally, the cutoffs are L/Q = 5 and ω Λ = √ L 2 + m 2 , while N f = 2N z = 64. The initial anisotropy was moderate, controlled by the parameters α = 2 and β = 4.
In Figure 6, we show the time evolution of τ n and τ in the unapproximated and in the classical schemes. τ n should be strictly constant 14 in both cases, since the conservation of particle number is not affected by the classical approximation (thus, this quantity is just used to monitor how well this conservation law is satisfied in the numerical implementation). A small difference between the two schemes is visible in the energy density. Given the conservation equation (35), view on this value is to recall that the screening mass in a φ 4 scalar theory at temperature T is m 2 scr = g 2 T 2 /24, while in Yang-Mills theory with 3 colors it is m 2 scr,YM = g 2 YM T 2 . Thus, if the two theories were compared at equal screening masses, one would have g 2 = 24 g 2 YM . Alternatively, if we compare the two theories at the same shear viscosity [34,35,36] to entropy density ratio, the scalar and gauge couplings should be related by A coupling g 4 = 50 in the scalar theory would correspond to a very small strong coupling constant αs ∼ 0.023 (conversely, αs = 0.3 would correspond to choosing g 4 ∼ 10 4 in the scalar theory). Note that if g 4 = 50 is the coupling at the scale Q, then the Landau pole of the φ 4 theory is at the scale µ = Q exp(16π 3 /(3g 2 )) ≈ 1844 Q -sufficiently above Q to justify a perturbative treatment. 14 Our discretization of the momentum integrals ensures exact conservation equations only if the time derivatives are evaluated exactly. The numerical resolution of the Boltzmann equation therefore also introduces an error that depends on the timestep ∆τ and on the details of the scheme used for the time evolution. With our implementation, the quantity (τ − ∆τ )n(τ ) is conserved with machine precision, and the expected conservation law is exactly recovered in the limit ∆τ → 0. this also indicates that the two schemes lead to different longitudinal pressures. Since the unapproximated scheme leads to a faster decrease of the energy density than the classical scheme, it must have a larger longitudinal pressure. This will be discussed in greater detail later in this section.

Bose-Einstein condensation
The initial condition that we have chosen corresponds to a large overpopulation, since [n −3/4 ] τ0 1. If the system were not expanding, we would expect the formation of a Bose-Einstein condensate. Figure 7 shows the particle density in the zero mode in the unapproximated and classical schemes. The onset of Bose-Einstein condensation is nearly identical in the two schemes, and a moderate difference develops at later times, that reaches about 20% at Qτ ∼ 10. The rather small difference between the two schemes for this quantity can be understood from the fact that the evolution of the condensate is governed by the region of small momenta, where the particle distribution is very large. We also see here a trend already observed in the isotropic case in ref. [32]: the classical approximation leads to more condensation than the unapproximated collision term. In fact, when the ultraviolet cutoff is large compared to the physical momentum scales, most of the particles tend to aggregate in a condensate in the classical approximation.
We mention in passing that our algorithm is not particularly well suited 22 for a detailed study of the infrared region. The fixed energy spacing of our discretization lacks resolution in the IR, and our treatment of the mass as fixed, rather than a self-consistently determined thermal mass, is another limitation.
On the other hand, the infrared has the highest occupancies, so classical methods are most reliable there. Therefore, lattice classical field simulations are much better suited for studying the infrared region. In particular our method is too crude to reveal the interesting scaling regimes found for instance in Ref. [2]. For this reason we will concentrate on quantities which are controlled by the higher-energy excitations, such as the components of the pressure.

Pressure anisotropy
More important differences between the two schemes can be seen in the behavior of the longitudinal pressure. In Figure 8, we display the time evolution of the ratio P L /P T . The beginning of the evolution is similar in the two schemes, with a brief initial increase of this ratio due to scatterings. Rapidly, the expansion of the system takes over and makes the ratio decrease, but at a pace slower than free streaming (indicated by a band falling like τ −2 ). The two schemes start behaving differently around Qτ ∼ 2, with the classical approximation leading to a faster decrease of the ratio P L /P T , approximately like τ −2/3 . In contrast, the unapproximated collision term seems to lead to a constant ratio at large times.
In Figure 9, we display the quantity β eff ≡ −τ d ln(P L /P T )/dτ as a function of time. If we parameterize and if the exponent β is slowly varying, then β eff gives the instantaneous value of this exponent. This figure is to be compared with Figure 2, where several scenarios for the behavior of this exponent have been presented. On this plot, we see that this exponent behaves quite differently in the two schemes, the asymptotic exponent being close to 2/3 in the classical approximation, while it is zero with the unapproximated collision term. Moreover, the exponent 2/3 does not appear to play any particular role when one uses the full collision term, since β eff does not spend any time at this value in this case, despite the large occupation number in this simulation. Therefore, the classical attractor scenario represented by the red curve in Figure 2 is not realized for this combination of initial condition and coupling. This computation also indicates that the condition f 1, that was regarded as a criterion for classicality, should be used with caution. In this example, it does not guarantee that the classical approximation describes correctly the evolution of the system. This condition is imprecise because f is in fact a function of momentum, and f 1 may not be true over all the regions of phase-space that dominate the collision integral.
Note that if the ratio P L / is approximately constant, (as is the case in the full calculation at large times) then we have and the overpopulation measure behaves as follows: For an isotropic system, δ = 1/3 and n −3/4 is a constant, and the Bose-Einstein condensate would survive forever (in our kinetic approximation where inelastic processes are not included). If the system remains anisotropic at large times, we have δ < 1/3 and n −3/4 decreases. Therefore, one expects that, if a condensate forms, it has a finite lifetime because the overpopulation condition will not be satisfied beyond a certain time. The final outcome should therefore be the disappearance of the condensate. The beginning of this process is visible in Figure 7 in the case of the unapproximated collision term. We have also investigated the sensitivity of our algorithm to the ultraviolet cutoff L on the longitudinal momentum. This cutoff may indeed have an important influence on the result since the typical p z of the particles in the system evolves with time due to the expansion. In Figure 10, we compare the results for L/Q = 5 and L/Q = 7. We observe that the difference between the two cutoffs is essentially the one inherited from the initial condition, i.e. the fact that the tail of the Gaussian in eq. (62) extends further when we increase the cutoff. The qualitative differences between the classical and full results are independent of the value chosen for this cutoff.
In Figure 11, we vary the coupling constant in order to see how the asymptotic behavior of the full solution is affected by the strength of the interactions. For the three values of the coupling, the ratio P L /P T reaches a minimum, whose value increases with the coupling. For the largest of the couplings we have considered (g 4 = 200, i.e. g ≈ 3.76), this ratio even shows a slight tendency to increase after a time of order Qτ ≈ 12.

More results using the DSMC algorithm
The deterministic algorithm we have used so far provides a direct resolution of the Boltzmann equation, but requires at each timestep the very time consuming computation of the collision integral. Moreover, it has a rather unfavorable scaling with the number of lattice points used in order to discretize momentum space. For this reason, we have also implemented a stochastic algorithm, the "direct simulation Monte-Carlo" (DSMC, described in the appendix E), where the distribution f is replaced by a large ensemble of simulated particles. By construction, energy, momentum and particle number are exactly conserved with this algorithm (provided the kinematics of the collisions is treated exactly). Its sources of errors are the limited statistics, and the reconstruction of the particle distribution from the simulated particles (this step requires a discretization of momentum space, which leads to some additional errors).
Before showing more results using this algorithm, we have first used it with the same initial condition already used with the deterministic method, in order to compare the two approaches. The outcome of this comparison is shown in Figure 12, and indicates a good agreement between the two methods. The differences, in the 10 % to 20 % range, can be attributed to the fact that the deterministic method simply disregards the particles with momenta higher than the lattice cutoffs. In contrast, in the DSMC method, there is no limit on the momenta of the simulated particles, and a discretization of momentum space is only used when reconstructing the distribution f from the ensemble of particles.
We have then used the DSMC algorithm in order to study the time evolution of the pressure ratio P L /P T in two situations:  Figure 12: Comparison of the deterministic method ("Direct") and the DSMC algorithm (see the appendix E), for the initial condition (62) with α = 2, β = 4 and g 4 = 50. In the DSMC case, the band is an estimate of the systematic error based on the different values of the pressure one obtains by including or not the particles from the "condensate" (since the definition of the condensate in the DSMC includes all particles in a small volume around p = 0).
First, we fix the value of the coupling constant at g 4 = 50, and we vary the prefactor f 0 in the initial condition given by eq. (62). The parameters α and β that control the initial anisotropy of this distribution are also held fixed, as well as the initial time Qτ 0 = 1. Although one may naively expect the agreement between the full and classical results to improve when f 0 is increased, Figure  13 shows that this is not the case. For the three values of f 0 considered in this computation, the classical approximation departs from the full result at roughly the same early time (or even a little earlier for the largest f 0 ). No matter how large f 0 , the ratio P L /P T becomes roughly constant at late times -or even slightly increases-in the full calculation, and decreases like (Qτ ) −2/3 in the classical approximation. Next, in a second series of computations, we have varied simultaneously the coupling constant and the initial occupation number in such a way that g 2 f 0 remains constant, g 2 f 0 ≈ 700. The parameters α and β controlling the initial anisotropy are the same as before. The resulting evolution of the ratio P L /P T is shown in the figure 14. In the classical approximation, the pressure ratio falls like (Qτ ) −2/3 , as already observed earlier. Note that we have represented only one classical curve, common to all the values of g 2 . Indeed, since the collision term in this approximation is homogeneous in f , one can factor out a prefactor f 2 0 and combine it with the g 4 of the squared matrix element, so that the classical dynamics is always the same if g 2 f 0 is held fixed. In contrast, the full dynamics does not posses this invariance, but appears to converge towards the classical result when g 2 → 0. At finite coupling, the agreement between the full and classical results is good only over a finite time window, that shrinks as g 2 increases. At moderate values of the coupling such as g 2 = 7 (corresponding to g 4 = 50, see Footnote 13), the unapproximated evolution departs from the classical one at a rather early point in time, and the exponent −2/3 does not play any particular role in the evolution of P L /P T . Two larger values of the coupling (g 2 = 45 and g 2 = 100) are also shown on this plot, but one should not take the corresponding results seriously. Indeed, scalar theory with such a large coupling is not really self-consistent, because the coupling runs very fast and the Landau pole is only a factor of 5 to 20 away.

Summary and conclusions
This paper started with the qualitative observation that large-angle out-of-plane scatterings are artificially suppressed by the classical approximation of the collision term in the Boltzmann equation with 2 → 2 scatterings, when the particle distribution is anisotropic, as is generically the case for a system subject to a fast longitudinal expansion. This kinematics is for instance realized in the early stages of heavy ion collisions. In order to quantify this effect, we have considered a longitudinally expanding system of real scalar fields with a φ 4 interaction in kinetic theory, and we have solved numerically the Boltzmann equation with elastic scatterings in two situations: (a) with the full collision term, and (b) in the classical approximation where one keeps only the terms that are cubic in the particle distribution.
This numerical resolution has been performed with two different algorithms. The first one is a direct deterministic algorithm, in which one discretizes momentum space on a lattice in order to compute the collision integrals by numerical quadratures (by assuming a residual rotation invariance around the p z axis, we could perform the azimuthal integrals analytically). The second method we have considered is a variant of the direct simulation Monte-Carlo (DSMC), in which the distribution is sampled by a large number of "simulated" particles.
The outcome of these computations is that the classical approximation is not always guaranteed to be good -even at a qualitative level-in situations where the occupation number is large. At moderate values of the coupling constant, the classical attractor scenario cartooned in Figure 2 is never observed, and the pressure ratio P L /P T becomes constant at late times or even increases without showing any sign of a τ −2/3 behavior in the full calculation (while a τ −2/3 behavior is indeed seen at late times in the classical approximation). Increasing the occupation number at fixed coupling does not make the classical approximation any better. The classical behavior, with τ −2/3 behavior in the longitudinal pressure, does emerge when one increases f 0 and decreases the coupling g 2 , keeping g 2 f 0 constant. But even in this case, the full quantum behavior deviated from the classical one when the occupancy at p ∼ Q, p z = 0 was still large.
This study indicates that the conventional criterion for classicality, f 1, is too simple in situations where f has a strong momentum dependence. If this condition is meant to be understood as f (p) 1 for all p, then it is not useful, because it is never realized. If instead one understands "f " as the maximal value of f (p), then this condition is necessary for the classical approximation, but by no means sufficient. Given this, the outcome of computations done in this approximation should be considered with caution unless confirmed by other computations performed in a framework that goes beyond this classical approximation. In addition, extrapolations of classical calculations from very weak coupling (where the classical and full calculations agree over some extended time window) to larger couplings must be taken with care, since the classical attractor behavior completely disappears at couplings that are still relatively small.
The limitation we have discussed in this paper is due to the missing f 2 terms when one uses the classical approximation in the collision term of the Boltzmann equation. Therefore, it does not affect the variant of the classical approximation mentioned in the last paragraph of Section 2, since the replacement f → f + 1 2 that one performs in this variant restores the f 2 terms. However, this variant suffers from a potentially severe sensitivity to the ultraviolet cutoff. For a nonexpanding system, one can mitigate this problem by choosing the cutoff a few times above the physical scale. In the appendix A, we show that this variant of the classical approximation describes isotropization in a non-expanding system much better than the plain classical approximation that has only the f 3 terms, without being much affected by the ultraviolet cutoff. The cutoff dependence of this variant of the classical approximation is much harder to keep under control in the expanding case, because the physical scales (and possibly the cutoff itself, depending on the details of the implementation) are time dependent. It also has much more severe problems when the collision term includes number-changing processes.
Let us end by mentioning the very recent work presented in ref. [37] where a similar study has been performed in the case of Yang-Mills theory, by applying the effective kinetic theory of ref. [38] to the study of a longitudinally expanding system of gluons. In this work, the authors also observe a rapid separation between the full evolution and the classical approximation, already for small couplings such as g 2 N c = 0.5.

A Anisotropic system in a fixed volume
In this appendix, we consider an anisotropic system of scalar particles in a fixed volume, in order to study how the classical approximation affects its isotropization. In this case, assuming that the particle distribution depends only on time but not on position, the left-hand side of eq. (5) reduces to so that the Boltzmann equation now reads with the collision term as given in eq (22). The evaluation of the collision term is identical to the case of a longitudinally expanding system, while the free streaming part of the equation is now completely trivial.
To obtain the numerical results presented in this appendix, we use again an initial condition of the form: with f 0 = 100, α = 2 and β = 1.2. Since the system is not expanding, the value of the initial time is irrelevant. We have taken Qt 0 = 0.1. The coupling constant is g 4 = 50 and the mass of the particles is set to m/Q = 0.1. The number of lattice spacings are set to N f = 2N z = 64, while the maximal values of p z and ω p are given by L/Q = 3 and ω Λ = √ L 2 + m 2 . In the case of an homogeneous non-expanding system, the conservation laws take a very simple form: that we use to monitor the accuracy of the numerical solution. Starting from the same initial condition given in eq. (69), we have studied the time evolution of the system with the full collision term (i.e. with both the cubic and quadratic terms in f ), in the classical approximation (with the cubic terms only), and in the "classical statistical approximation" (CSA) that amounts to including the zero point vacuum fluctuations in the corresponding classical field approximation. See Section 2.6 for more details on the different classical schemes.
In Figure 15, we compare the time evolution of the particle density in the condensate, n c , for the three schemes. It appears that the three schemes agree qualitatively, and even semi-quantitatively, for the evolution of this quantity. The onset of condensation is almost exactly identical in all the schemes, while the final values of n c differ by about 10%. In agreement with the observations of ref. [32], the classical approximation tends to overpredict the fraction of condensed particles, while the classical statistical approximation tends to underpredict it. Next, we consider the time evolution of the ratio between the transverse and longitudinal pressures. In Figure 16, we plot the time evolution of P T /P L − 1 in the three schemes. Starting from a nonzero value dictated by the momentum anisotropy of the initial condition, this quantity is expected to return to zero as the particle distribution isotropizes. After a short initial stage during which the three schemes are nearly undistinguishable, we observe that the unapproximated scheme and the classical statistical approximation lead to almost identical time evolutions for this quantity, while the classical approximation isotropizes at a much slower pace. This is consistent with the argument exposed in the introduction, according to which the terms quadratic in f are essential for out-of-plane scatterings in an anisotropic system.

Expression as an elliptic integral
In eq (20), we have encountered an integral involving the product of four Bessel J 0 functions, which, to the best of our knowledge, does not seem to be known in the literature.
In this appendix, we derive an exact expression of this integral in terms of an elliptic function K, starting from the known expressions of similar integrals with two J 0 Bessel functions (see ref. [39]-6.512.8): and three J 0 Bessel functions (see ref. [40]): A(p 1 , p 2 , p 3 ) is the area of the triangle whose edges have lengths p 1 , p 2 and p 3 (if such a triangle does not exist, then the integral is zero). We can therefore recast eq. (71) into the following expression: .
Recall that the area of a triangle in terms of the lengths of its edges is given by the following formula, Thus, eq. (71) can be expressed as where we denote α ij ≡ (p i + p j ) 2 and β ij ≡ (p i − p j ) 2 . The new integration variable is x ≡ t 2 . The range of integration on x is restricted by the fact that the arguments of the two square roots must both be positive: Since it involves only a square root of a fourth degree polynomial, the integral in eq. (76) is an elliptic integral, which can be reduced to a combination of Legendre's elliptic functions.

B.2 Expression in terms of the elliptic K function
The boundaries of the integration range in eq. (76) are two of the roots of the polynomial Let us call r 1 and r 2 these two roots, respectively the lower and upper bound. And for definiteness, let us call r 3 and r 4 the remaining two roots, arranged so that r 3 < r 1 < r 2 < r 4 , i.e. r 3 ≡ min(β 12 , β 34 ) , r 1 ≡ max(β 12 , β 34 ) , Eq. (71) is thus equivalent to Let us now perform the following change of variables: x = r 1 cos 2 θ + r 2 sin 2 θ .
The above integral becomes with α 2 ≡ (r 4 − r 2 )/(r 2 − r 1 ) > 0 and β 2 ≡ (r 1 − r 3 )/(r 2 − r 1 ) > 0. This integral can be expressed in terms of the complete Legendre elliptic function of the first kind (defined for z ∈ [0, 1)), which leads to the following compact expression 15 The function K(z) can be evaluated efficiently with a simple algorithm based on the arithmetic-geometric mean: (85) Figure 17 shows a comparison of a direct numerical evaluation of the integral in eq. (71) with the formula (84). Note that this quantity becomes singular for special configurations of the p ⊥i 's (in particular, when their values allow the vectors to become collinear). Near these points, the direct method is inefficient because of the very slow convergence of the integral. In contrast, eq. (84) is much better because the algorithm described in eqs. (85) converges to a very accurate result in only a few iterations 16 . Moreover, this method does not require the evaluation of any transcendental function.

C Integration domain for the collision term
When we discretize the collision integral of eq. (22), the domain of integration for the energy variables ω p3 and ω p4 is the one represented in Figure 18. When 15 From this formula, one can check the identity I 4 (p 1 , p 2 , p 3 , 0) = (2πA(p 1 , p 2 , p 3 )) −1 , that one expects from eqs. (71) and (73). 16 The number of accurate digits doubles at every iteration.    Figure 19: Discretization of the integrals over p z3 and p z4 , with the values of the quadrature weights of each point. i 1 = 1 or i 1 = N f , this domain becomes a triangle. In this case, the 6 points represented in green merge pairwise to form the summits of the triangle. The quadrature weight at these points becomes 1 2 × 1 2 × 1 2 = 1 8 . But we do not need to handle this case by hand since the formula (50) gives the correct weights in all cases. Likewise, we have represented the integration domain on the longitudinal momenta p z3 and p z4 in Figure 19, with the quadrature weight of each point. Since we have assumed that p z1 > 0, the index j 1 is positive.
We need also to specify the integration domains for the terms C 1c↔34 p1 [f ] and C 12↔c4

p1
[f ] that describe collisions between a particle from the condensate and a particle at non-zero momentum, whose expressions are given in eq. (29). When p z1 > 0, the integration over ω p4 in C 1c↔34

p1
[f ] is discretized as a sum over the following discrete points: while for the longitudinal momentum we have: For C 12↔c4 p1 [f ] given in eq. (29), the integration domain for the energy ω p4 is 38 while for the longitudinal momentum p z4 we have Finally, we need also to specify the discrete domains in the right hand side of the equation (31) for the evolution of the particle density in the condensate. The sum over the energies ω p3 and ω p4 is over the following set of points: In this section we derive the weights (57) that enter into eq. (52). This derivation can be done in the case where there is no Bose-Einstein condensate, since the particles in the BEC carry zero momentum and are therefore not affected by free streaming. Let us start with the particle density defined in eq. (53). From the conservation of the number of particles, it should obey τ ∂ τ n nc = −n nc .
Let us recall that eq. (52) is only valid for j > 0. Therefore, we first need to rewrite the particle density (53) by using the parity in p z (i.e. h[i, −j] = h[i, j]) and the fact that w z [−j] = w z [j] As one can see on Figure 5, particles can hop to the j = 0 line from the left or from the right. Therefore, for j = 0, eq. (52) can be rewritten as follows: Then , by summing eqs. (52) and (92) on the indices i, j, we obtain (the indices of the last two terms have been shifted) The right-hand side must therefore be equal to that of eq. (91), for every h[i, j]. This leads first to and Next, using the definition for the total energy (54) and the longitudinal pressure (55), we can use Bjorken's law vi. Increment the time t → t + ∆t and return to the step i.
In order to illustrate this method and its difference with the deterministic method used in the rest of this paper, we have applied it to the case of a spatially homogeneous non-expanding system. In this case the Dirac function is regularized by δ(p) → of the condensate in the DSMC method also includes the particles in a small sphere of radius p min around p = 0, the transition of condensation appears less sharp than in the deterministic method where the density n c is defined as the coefficient of the delta function δ(p). By decreasing the value of the p min used in this definition, the transition in the DSMC simulation appears to become sharper, to finally become very close to the transition observed in the direct method. Alternatively, this interpretation can be further checked by adding to the n c of the deterministic method the contribution of the non-condensed particles contained in this small sphere.