The Boltzmann equation in classical Yang–Mills theory

We give a detailed derivation of the Boltzmann equation, and in particular its collision integral, in classical field theory. We first carry this out in a scalar theory with both cubic and quartic interactions and subsequently in a Yang–Mills theory. Our method does not rely on a doubling of the fields, rather it is based on a diagrammatic approach representing the classical solution to the problem.


Introduction, motivation and the Boltzmann equation
Transport phenomena in QCD matter have been the subject of extensive research over the last three decades. Particular attention has been paid to calculating quantities like conductivity, viscosity and baryon diffusion [1][2][3][4] or the relaxation of colorful excitations [5][6][7][8][9][10] in a weakly coupled Quark-Gluon Plasma (QGP). A key element in such studies has been the use of kinetic equations which are of the Boltzmann type. The Boltzmann equation is an equation which describes the time evolution of occupation numbers. An occupation number is a dimensionless quantity defined as the number of particles of a given species per unit phase space and divided by the number of choices for each possible discrete degree of freedom. For example, in a SU (N c ) pure gauge theory one divides by 2(N 2 c − 1) for the polarizations and colors of the gauge bosons to which we shall refer as gluons. The Boltzmann equation for the gluon occupation number f ( p, x, t) reads where we have used the compact notation f p = f ( p, x, t) since the integrand is local in both x and t and defined in general the integration measure Energy-momentum conservation in Eq. (1.1) is explicit, while the scattering amplitude squared |M| 2 YM for the process p 2 p 3 → pp 1 is summed over initial and final colors and polarizations and is given below in Eq. (3.19). Each of the two terms in the square bracket in Eq. (1.2) has an intuitive interpretation. The first is a gain term proportional to f p 2 f p 3 , with p 2 and p 3 disappearing to create p and p 1 , while 1 + f p 1 and 1 + f p are Bose enhancement factors. Similarly, the second is a loss term describing the disappearance of p and p 1 in order to create p 2 and p 3 . Notice also that this square bracket vanishes when occupation numbers are given by the Bose-Einstein distribution. Further aspects of this collision integral will be discussed in the next sections.
A valid question that one immediately asks is how such a kinetic equation can be derived from first principles, i.e. from the underlying quantum field theory. Indeed, this was first addressed long time ago in non-relativistic quantum field theory [11]. Using the Schwinger-Keldysh formalism and writing Dyson-Schwinger equations for the propagators, an appropriate truncation supplemented with a gradient expansion led to the non-relativistic version of the Boltzmann equation given above. Notice that in such a limit the Bose enhancement factors in the collision integral are absent and the collision integral vanishes when occupation numbers are given by the Maxwell-Boltzmann distribution. Using similar Green's function techniques in relativistic quantum field theories, the Boltzmann equation was derived in [12] for scalar fields, in [13] for charged scalar fields and in [14] for nuclear matter described by the Walecka model. A somewhat different derivation based on resuming ladder diagrams, again in a scalar field theory, was given in [15], while kinetic equations for colorful excitations in a weakly coupled QGP were obtained in [9,16,17] by performing gauge covariant gradient expansions. For both a pedagogical introduction and an overview we refer the reader to [10,18].
Typically, the essential assumptions for arriving at such a kinetic equation are two. First one needs that occupation numbers do not become very large; for example in QCD one needs f p 1/α s while in a scalar theory with quartic interactions (λφ 4 theory) this constraint would be f p 1/λ. This is necessary, since otherwise a description using onshell scattering of individual particles no longer makes sense as the time between scatterings is too short for an on-shell approximation to be valid. Second one has to assume that there are no large wavelength modes comparable to the mean free path, otherwise one has to treat them in a suitable way.
Here we would like to study the conditions under which bulk matter can be described by a Boltzmann equation with a collision term given by elastic scattering, but also under the additional assumption that the physical system is classical 1 . Then the extra condition f p 1 is required in order to have the possibility of a quantum-classical correspondence, but when the coupling is sufficiently small there is a parametrically large window in which a kinetic description via a Boltzmann equation should be valid. In fact such an observation and the corresponding derivation have been already done a few years ago in the context of a λφ 4 theory [21] (see also [22]). In that work, the starting point of the analysis was a doubling of the fields, a method which has been naturally used for the corresponding quantum problem where separate fields are needed for time evolution in the direct amplitude and the complex conjugate amplitude. However, when occupation numbers are large one combination of the fields, π in [21], becomes a variable of constraint and the functional integration over π requires the other independent combination 1 A different connection between the classical approximation to statistical field theory and the transport theory appears in studies of baryon number violation via topological transitions in hot QCD; in that context, the quantum Boltzmann equation for the relaxation of colorful excitations has been used to construct a classical effective theory for the "ultrasoft" modes responsible for the topological transitions [7,17,19,20]. of fields, φ, to obey the classical equations of motion of the λφ 4 theory. Thus, although there is only one dynamical variable in the discussion given in [21], the constraint variable appears explicitly in the perturbative classical calculation of the Boltzmann collision term.
There are two major differences between the current work and the one in [21]. The first is that we simply use a different method which does not rely on the doubling of the fields; we solve classical equations of motion, with retarded boundary conditions as appropriate to the problem, in which only one field evolves and interacts. Occupation numbers are not defined in terms of Green's function, as usually done in the quantum analyses and in that of [21]. Instead we start from the "canonical" definition that f p should be proportional to a * p a p where a * p and a p are the classical analogs of creation and annihilation operators, i.e. the coefficients in the expansion of the classical field in plane waves. In this language it is clear how the constraint f p 1 emerges, since in the classical treatment we consider these expansion coefficients as numbers and not as operators, thus effectively ignoring all possible commutators. Now we can follow the classical time evolution of the field coefficients and, in turn, that of the occupation numbers.
The second difference with respect to [21] is that we extend the analysis to the case of a Yang-Mills theory. In order to efficiently deal with the latter, we shall first consider a scalar theory with both cubic, gφ 3 , and quartic, λφ 4 , interactions. Then the study of the Yang-Mills theory becomes much easier since the topology in the diagrammatic expansion is the same, with the only additional complication being the introduction of spin and color degrees of freedom. Our calculations, using classical field equations as already stressed, are given as the first terms in a power series in g 2 and λ in the scalar theory and in g 2 in the Yang-Mills theory. They agree with the corresponding quantum field theory result so long as occupation numbers satisfy f p 1 and after ensemble averages (whose particular details should not matter when the constraints in the occupation numbers are satisfied) over the initial conditions are performed in both the classical and quantum approaches. Thus, we shall eventually arrive at the collision integral in Eq. (1.2), but it will contain only the cubic in f terms and not the quadratic ones, cf. Eq. (3.25). The equilibrium limit in that equation is now given by f p = kT /E p , which is clearly the large occupation limit of the Bose-Einstein distribution occurring when E p kT . In order to make our discussion as simple as possible we have made a number of assumptions: (i) We suppose that the elements of our initial ensemble of field configurations are homogeneous in space. This assumption is not really necessary, but it simplifies our task considerably. What one must actually assume is that inhomogeneities occur on a scale large compared to the wavelengths dominating the problem and this is sufficient to get an effective momentum conservation, e.g. the δ (3) ( p) emerging in Eq. (2.17). When such spatial inhomogeneities are present they trivially give rise to the drift term v p · ∂ f p /∂ x which appears in the Boltzmann equation in Eq. (1.1) and combines with ∂ f p /∂t term to form the natural "convective" derivative. (ii) We assume the absence of long range coherent fields which would give rise to the term F ext · ∂ f p /∂ p in Eq. (1.1). (iii) We finally suppose that our initial fields ensemble does not have long range coherences in wavelengths so that Eq. (2.4), which defines the occupation numbers, is appropriate. Similar assumptions were made in the analysis of [21], however, other possibilities are available as we now discuss.
The above assumptions are generally satisfied in recent studies of scalar field theories and their simulations [23][24][25]. However, in simulations of Yang-Mills theories this is not always the case. On the one hand, in [26,27] the initial conditions are very much as we have taken them and one expects that after a short time, allowing occupation numbers to become smaller than 1/α s , the classical field theory simulations should agree with the Boltzmann equation. Indeed, this seems to be the case, as the results in [26,27] are very close to the Boltzmann-based description given in [28]. On the other hand, the recent simulations in [29] begin with long range coherent fields and thus Eq. (2.4) is not satisfied. At this point it is not clear at what time the classical field evolution of [29] would admit an equivalent description via a Boltzmann equation.
In Sect. 2 we do the derivation for the scalar theory with gφ 3 and λφ 4 interactions. The calculation is based on suitable Feynman rules which allow for a diagrammatic solution of the classical equations of motion. We have separated the calculation in three subsections in which we calculate in great detail the λ 2 , the λg 2 , and the g 4 terms respectively. Each of the aforementioned terms contains all the gain and loss terms of the collision integral. Then, in Sect. 3, we give the derivation for a Yang-Mills theory by paying special attention to the points that require extra treatment compared to the scalar theory case.

Scalar field theory with cubic and quartic vertices
Let us start by considering a massless scalar field theory with cubic and quartic interactions in D = 4 dimensions. The action is given by and while the coupling λ is dimensionless, the coupling g has mass dimension 1. In this work, and in view of the perturba-tion theory to follow, we shall assume that λ and g 2 /M 2 are of the same order, where M is a typical mass scale for the scattering processes to be taken into account. In general, we can decompose the real classical field φ according to and where p is an on-shell four-momentum so that p · x = E p x 0 − p · x and E p = | p|. Since we have an interacting field theory, the coefficients a p and a * p are generally time-dependent. However, the Boltzmann equation is valid when the typical collision time is much smaller than the time between two collisions. Thus, even though we will assume that a p is time-dependent, we will take this dependence to be much slower than that of the plane wave in Eq. (2.2). This allows us to invert Eq. (2.2) and express a p in terms of the field φ as In the case of a homogeneous medium it is natural to define the occupation number f p , a dimensionless quantity, as with the shorthand notation δ and where the brackets stand for the ensemble average. We aim to find the time evolution of the occupation number in the classical theory and therefore we need to determine the corresponding evolution of the coefficients a p and the field φ. The classical equation of motion of φ clearly reads with the convention x = ∂ 2 0 − ∇ 2 x and where we have defined for our convenience the "current" J . Let us now split the full interacting field φ according to where φ (0) is the free field, i.e. it satisfies the homogeneous version Eq. (2.5), while δφ is the modification arising from the presence of interactions; it satisfies Eq. (2.5) and thus can be formally written as In the above is the free propagator of the scalar field and is determined by The solution to the above is where → 0 + so that the propagator is proportional to (x 0 ), as is straightforward to check by performing the integration over k 0 . More precisely, one finds Therefore, the propagator in Eq. (2.9) is the retarded (or causal) one, since this is the natural choice when initial conditions (that is, φ (0) ) are given. For later use let us note that this retarded propagator can also be written as where P stands for principal value and therefore one has a clear separation of the real and imaginary contributions to the propagator. Now, in analogy to Eq. (2.6) we can split the coefficient a p as a p = a (0) p + δa p (2.12) and using the form of the propagator given in Eq. (2.10) just above we easily find that the piece δa p generated by the interactions is given by The corresponding change in the occupation number reads where, in writing the first term on the r.h.s. of the above, we have anticipated that it will be proportional to δ (3) p p like the l.h.s. Finally, by taking a time derivative we arrive at where, with a slight notational abuse, δȧ p stands for the time derivative of δa p . We shall refer to the two terms in the r.h.s. of Eq. (2.15) as the crossed and diagonal terms, respectively. In general one cannot solve Eq. (2.7) and/or Eq. (2.13); that would be equivalent to solving the full nonlinear classical problem, which is in any case beyond our goals. What we shall do is to assume that the correction δa p is small compared to a (0) p and perform a calculation to first non-vanishing order in λ ∼ g 2 /M 2 . Eventually this translates to imposing the condition that occupations numbers do not get large, more precisely f p 1/λ. Recalling that the classical approximation to the problem also requires f p 1, we see that there is a parametrically large window of validity for the "classical" Boltzmann equation, so long as the couplings are sufficiently small.

2.1
The λ 2 terms and the Feynman rules for classical diagrams in the scalar theory To illustrate the procedure, we shall first do a step-by-step calculation for the λ 2 contribution to the diagonal term in Eq. (2.15), which simply means that we need to find the order λ contribution to δa p . Since the current in Eq. (2.5) is already of order λ we can substitute the full field φ with its free part φ (0) . Next, for reasons to become apparent in a while, let us consider the following particular term in where p 1 , p 2 , and p 3 are on-shell four-momenta and with the combinatorial factor 3 coming from the number of ways we can pick the required product of field coefficients out of [φ (0) ] 3 . Now we can integrate over y to get 3 . Notice that we have dropped the superscript (0) from the expansion coefficients, since this is allowed to the level of accuracy and in order to have a more economical notation. Furthermore, let us point out that at this stage energy is not conserved at the vertex. The y 0 time integration is unbounded for large negative values and we make it convergent via the "adiabatic" prescription E → E − i with → 0 + to find From the above "direct amplitude" (DA) it is straightforward to construct its time derivative δȧ p and the "complex conjugate amplitude" (CCA) δa * p . When forming δa * p δȧ p we encounter a six-point correlator of the field coefficients and since the system is dilute we will assume that it factorizes to a product of two-point functions, that is, to a product of occupation numbers. More precisely, we assume the ensemble average and since we integrate over all momenta one immediately sees that the two terms in the above will eventually contribute the same to the final result. Using the δ-functions arising from the ensemble average in Eq. (2.19) one can readily perform all the integrations over the primed momenta in the product δa * p δȧ p . Then the δ-function corresponding to momentum conservation in the CCA becomes δ (3) and after also using momentum conservation in the DA it finally gives a factor δ (3) p p as expected (cf. the discussion after Eq. (2.14)). Now E becomes the same in the DA and in the CCA and we have which is the required energy conservation. Now we put everything together in Eq. (2.15) to finally arrive at the λ 2 gain terṁ where p = p + p 1 − p 2 − p 3 with all four-momenta being on-shell and where we have adopted the compact notation introduced in Eq. (1.3) for the integration measure. Let us note here that it is only the choice made in Eq. (2.16) for the field coefficients which leads to energy conservation. Any other combination, e.g. an a * a * a term, will lead to complex exponentials with uncompensated energy differences. Such exponentials will average to zero at large times, since the time scales describing variations in the Boltzmann equation are supposed to be very large compared to the typical interaction times. λ 2 is simply the amplitude squared |M( p 2 p 3 ; pp 1 )| 2 in the λφ 4 theory and Eq. (2.21) acquires a natural interpretation as a gain term arising from a 2 → 2 scattering. The integrand is naturally proportional to the occupation numbers of the incoming particles f p 2 and f p 3 , while f p 1 appears as a Bose enhancement factor. The (square of the) Feynman diagram related to the term we have just calculated is shown in Fig. 1.
Let us now establish some Feynman rules for the classical problem at hand in order to systematize the calculation for the remaining terms. For any diagram in the DA we have the following momentum space rules: Fig. 1 The λ 2 contribution to δa p δa * p , cf. Eq. (2.18). A circled cross stands for an external insertion while the open line corresponds to the momentum measured. The ensemble average will set p 1 = p 1 , p 2 = p 2 and p 3 = p 3 (or p 2 = p 3 and p 3 = p 2 ), while momentum conservation in both the DA and in the CCA will lead to p = p • Assign a factor 1/ h p from the definition of δa p .
• Assign a factor −ig for each cubic vertex and a factor −iλ for each quartic one. • Divide by the symmetry factor. The maximum such factor we will come across is 2; this will take place when two field coefficients of the same type, that is, two a's or two a * 's, are connected to the same vertex. • Impose three-momentum conservation at each vertex.
• Assign an overall factor (2π) 3 δ (3) ( p) where p is the sum over all external three-momenta in which the momentum p and the momenta associated with a * 's are taken with a positive sign, while the momenta associated with a's are taken with a negative sign. • Impose energy conservation at all, but one (see next rule), vertices.
→ 0 + at the vertex which connects to the measured occupation factor. E is the energy imbalance at the vertex, and thus also that of the full diagram, with E p taken with a positive sign.
• Use the retarded propagator R (k) = i/(k 2 + i k 0 ), with → 0 + , for each internal line. The four-momentum k should flow towards the measured occupation factor. Equivalently, one can use the advanced propagator if the four-momentum k is taken to flow away from the measured occupation factor.
• Integrate according to d 3 p h p a * p or d 3 p h p a p for each external line, but not for the measured particle.
We stress that these rules are just a convenient representation of the perturbative solution to the classical problem. It is trivial to check that they lead to Eq. (2.18) when considering the DA in Fig. 1.
Next, we shall use these Feynman rules to calculate the remaining λ 2 terms. These come from the crossed term in Eq. (2.15) and it is clear that now we need to compute δa p to order λ 2 . To this order, the two diagrams which will eventually satisfy energy conservation are shown in Fig. 2. As we shall see, the Fig. 2.a leads to the loss terms in the Boltzmann equation while 2.b leads to a gain term. Even though it is not necessary, let us mention, just for illustrative purposes, that such diagrams arise from the current J (y) expanded to order λ 2 which can easily be found to be where we have dropped the superscript (0) in the field φ. Now one would need to expand all the free fields in plane waves as before, but as explained above it is more convenient and much less tedious to directly use the Feynman rules. We readily see that the Fig. 2.a gives 6 , and k = p 5 + p 6 − p 4 . The symmetry factor 2 in the denominator comes about because the diagram remains invariant under the exchange of the legs corresponding to momenta p 5 and p 6 . Differentiation w.r.t x 0 cancels the energy denominator and multiplication with a * p (cf. Eq. (2.15)) leads again to a product of six field coefficients. As in Eq. (2.19) we assume that the six-point correlator factorizes into a product of occupation numbers, that is, The factor of 2 comes because p 5 has to be contracted with either p 1 or p (and, correspondingly, p 6 with either p or p 1 ) and both terms contribute equally. The δ-function in the integrand of Eq. (2.23) reduces to δ (3) p p , and then E vanishes and k becomes p + p 1 − p 2 . Furthermore, making use of Eq. (2.11) we have which expresses energy conservation. Notice that due to the three δ-functions in Eq. (2.24), there are only two three-momentum integrations to be done which means the δ-function of the three-momentum conservation has been already implicitly used. To comply with the notation of Eq. (2.21) one can re-insert an integration over the momentum k, which we rename to p 3 , accompanied by δ (3) Then by putting everything together in Eq. (2.15) we arrive at the order λ 2 loss termṡ where, as in Eq. (2.21), p = p + p 1 − p 2 − p 3 with all four-momenta on-shell. Notice that we have been allowed to  6 , and k = p 4 + p 5 − p 6 . The symmetry factor 4 in the denominator comes about because the diagram remains invariant under the exchange of the legs corresponding to momenta p 1 and p 2 and the exchange of the legs corresponding to p 4 and p 5 . The six-point correlator factorizes into a product of occupation numbers according to a p 1 a p 2 a * p 4 a * p 5 a p 6 a * p → 2δ (3) where the factor of 2 arises because one can set p 4 = p 1 , p 5 = p 2 or p 4 = p 2 , p 5 = p 1 . The momentum k becomes p 1 + p 2 − p. For the propagator, which is advanced since we took the momentum to flow away from the measured occupation factor, we have This is the point where the two diagrams in Fig. 2 differ from each other. Compared to Eq. (2.25) the sign in Eq. (2.29) has changed and therefore Fig. 2.b leads to a gain term. As before we insert an integration over the momentum k, which we rename to p 3 , accompanied by δ (3) ( p+ p 3 − p 1 − p 2 ) and we immediately let p 1 ↔ p 3 . We put everything together in Eq. (2.15) to arrive at the order λ 2 second gain terṁ where p is as in Eqs. (2.21) and (2.26). Again, as already explained below Eq. (2.26) for the corresponding loss term, the Fig. 2.b eventually acquires an interpretation as 2 → 2 scattering. The integrand is proportional to the scattering amplitude squared λ 2 and to the occupation numbers f p 2 and f p 3 of the "incoming" momenta while f p is a Bose enhancement factor. Adding all the λ 2 contributions from Eqs. (2.21), (2.26), and (2.30) we arrive in fact at the Boltzmann equation in the classical φ 4 theory, that is, (2.31)

The g 4 terms
Let us turn our attention to contributions arising solely from the cubic vertices, i.e. the g 4 terms. What is non-trivial, compared to the λ 2 terms, is that now the amplitude squared |M| 2 depends on the kinematics. This dependence, containing the well-known s, t, and u diagrams, should come out from our calculation. Before writing down the diagrams, and focusing first on the diagonal term δa * p δȧ p in Eq. (2.15), we give for completeness the current J (y) to order g 2 ; a single iteration leads to Fig. 3 The g 2 contributions to δa p leading to a gain term in the Boltzmann equation The Feynman diagrams for δa p , which in the end will contribute to the Boltzmann equation, are shown in Fig. 3. In analogy to the corresponding λ term [cf. Eq. (2.18)] we need a product of the type a * aa, and since a * can originate either from φ(y) or from φ(z) we have the two distinct diagrams in Fig. 3. Using the Feynman rules we can combine both diagrams into with E and p as in Eq. (2.18). The two diagrams differ only in the symmetry factors (1/2 and 1 respectively) and in the argument of the retarded propagator. Eq. (2.33) is very similar to Eq. (2.18) with the only difference being the presence of a propagator for each of the two terms. In fact, the only role of these propagators is to lead to the proper form of |M| 2 in the gφ 3 theory. Therefore the calculation is almost identical to the one following Eq. (2.18). In particular, notice that the real part of the propagators, since they are in general off-shell, does not play any role in the computation of the diagrams under current consideration and the energy conservation will emerge as in Eq. (2.20). We just need to be careful to pick-up the correct arguments of the propagators after the contractions between the DA and the CCA due to the ensemble average. Defining the Mandelstam variables it is just a matter of simple bookkeeping to find the propagator products after taking the ensemble average of δa * p δȧ p . For p 1 = p 1 , p 2 = p 2 , and p 3 = p 3 (with the prime denoting momenta in the CCA) we have Putting everything together and noticing that one can let 2/t 2 → 2/u 2 and 2/st → 2/su inside the integrand we find the gain terṁ Considering now the crossed term a * p δȧ p in Eq. (2.15), one needs to calculate δa p to order g 4 . After straightforward iterations one finds that the current J (y) to this order reads (in a compact notation where repeated coordinates are integrated over) In Fig. 4 we show the six diagrams contributing to δa p . All corresponding expressions are very similar to Eq. (2.23) with the extra element of having two more propagators. We have (2.39) with E and p as in Eq. (2.23).
At this point it is appropriate to say that only a propagator with argument the sum of three external momenta will have a real part leading to conservation of energy. In fact, we have already used this property when considering the propagators in Eq. (2.33); none of the two propagators there acquired a real part. Moreover, this is also the reason that no diagram coming from the last term of the current in Eq. (2.38) contributes to the Boltzmann equation; any propagator in such a diagram will have as an argument the sum of either two or four external momenta, as one can easily verify by simply drawing it. Thus, energy conservation will emerge out of Eq. (2.39), as in Eq. (2.25), from the propagator R ( p 5 + p 6 − p 4 ) and the only extra work we have to do is to carefully calculate the arguments of the remaining propagators in the square bracket in Eq. (2.39) after taking the ensemble average of a * p δȧ p and without worrying about their real parts. One has to always identify p 4 with p 2 , while there is the possibility to choose p 5 = p 1 , p 6 = p or p 5 = p, p 6 = p 1 . It is an easy exercise to verify that the sum of propagator products in the square bracket in Eq. (2.39) becomes R R → −2 (2.40) Following now the exact same steps as in the case of the corresponding λ 2 term, and noticing in particular that the integrand is still invariant under p 2 ↔ p 3 , and thus under t ↔ u, so that we can let 2 3 , we arrive at the g 4 loss termṡ (2.41) Regarding the second g 4 gain term, we can draw four diagrams with the external lines a p 1 a p 2 a * p 4 a * p 5 a p 6 . Then, in analogy to the computation performed for Fig. 2.b, it is not hard to convince ourselves that we get a contributioṅ Now we put together all the g 4 contributions from Eqs. (2.37), (2.41), and (2.42) to arrive at the Boltzmann equation in the classical φ 3 theory, that is, The λg 2 terms and the Boltzmann equation for the full scalar theory Finally, in order to complete the derivation of the Boltzmann equation in the full scalar theory, i.e. with both cubic and quartic vertices, we need to compute the terms of order λg 2 . The first gain term emerging from the product δa * p δȧ p is rather easy to obtain since we already have the λ and g 2 contributions to δa p as given in Eqs. (2.18) and (2.33), respectively. Compared to the corresponding calculation of the λ 2 and g 4 terms, the only difference in this mixed term is coming again from the propagators which after the ensemble average give (2.44) In the above we have user for one more time our freedom to let 1/t → 1/u due to the invariance of the integrand in the subsequent integrations. We finally find the gain terṁ Regarding the crossed term a * p δȧ p term in Eq. (2.15), we have to compute δa p to order λg 2 . After straightforward iterations we find that the current J (y) to this order is given by In Fig. 5 we present the five diagrams contributing to δa p to order λg 2 . All corresponding expressions have similar structure to that of Eqs. (2.23) and (2.39), more precisely we have with E and p as in Eq. (2.23). Energy conservation will come from the propagator R ( p 5 + p 6 − p 4 ) as in the respective λ 2 and g 4 terms. Taking the ensemble average in the product a * p δȧ p we will identify p 4 with p 2 , while there is the possibility to choose p 5 = p 1 , p 6 = p or p 5 = p, p 6  Now we copy the same steps as in the case of the corresponding λ 2 and g 4 terms to arrive at the λg 2 loss terms, (2.49) Concerning the second λg 2 gain term, we can draw four diagrams with the external lines a p 1 a p 2 a * p 4 a * p 5 a p 6 . Then, in analogy to the previous respective computations, we get the expected contributioṅ It is trivial to add Eqs. (2.45), (2.49), and Eq. (2.50) to get the total λg 2 contribution. By furthermore adding the total λ 2 and g 4 expressions given in Eqs. (2.31) and (2.43), we come to the Boltzmann equation for the full scalar theorẏ where we have defined the scattering amplitude squared of the full scalar theory (2.52) Here we would like to stress that the specific combination of the occupation numbers in Eq. (2.51) and the scattering amplitude squared of the scalar theory have emerged as a result of our calculation. Let us also notice that a factor 1/2 in front of the integral in Eq. (2.51) is a symmetry factor due to the fact that particles 2 and 3, whose momenta are integrated over, are identical.
Furthermore, notice that the explicit form of |M| 2 φ as given in Eq. (2.52) was derived in detail in the context of this scalar field theory. In the Yang-Mills case, which follows in the next section, we shall not derive the respective amplitude squared |M| 2 YM , since this is a standard, albeit not trivial, textbook calculation. However, we shall of course show that |M| 2 YM emerges in all terms in the Boltzmann equation and this is sufficient for our proof. Thus, it is useful to reflect back and see how we arrived at |M| 2 φ in this section. This is straightforward for the diagonal gain term; combining Eqs. (2.18) and (2.33) we see that M φ ( p 2 p 3 ; pp 1 ) appears in the integrand in the DA. Similarly M * φ ( p 2 p 3 ; pp 1 ) appears in the CCA and after squaring and performing the ensemble average we arrive at |M( p 2 p 3 ; pp 1 )| 2 φ . Regarding the crossed term it is enough to look, for example, in the loss terms and a first discussion has already appeared below Eq. (2.26) in the λφ 4 case. Putting together Eqs. (2.23), (2.39), and (2.47) we see that M φ ( p 2 k; pp 1 )M * φ ( p 4 k; p 5 p 6 ) appears in the DA. After multiplying with the CCA taking the ensemble average and using the fact that k is put on-shell according to Eq. (2. 25) we arrive again at |M( p 2 p 3 ; pp 1 )| 2 φ (cf. the renaming of the momentum k below Eq. (2.25)).

Yang-Mills theory
Now we would like to extend our analysis to the Yang-Mills theory in D = 4 dimensions. Even though we will keep the number of colors N c arbitrary, we shall refer to the gauge bosons as gluons. The topology of the diagrams is the same as that in the full scalar theory studied in Sect. 2 and the extra complications come only from the color and spin structure of the diagrams. The Yang-Mills action in an axial gauge reads with the field strength and where f abc are the familiar structure constants of the SU (N c ) group. In general, n μ and ξ are arbitrary in Eq. (3.1), but for our convenience we shall consider the light-cone gauge defined by the conditions n μ n μ = 0 and ξ → 0. Introducing the polarization vectors ε λ μ ( p) for the two transverse (and physical) gluon polarizations, and which satisfy p · ε λ ( p) = n · ε λ ( p) = 0, we can expand the gauge field as Assuming a λa p is slowly varying and using the orthogonality property of the polarization vectors, i.e. ε λ ( p) · ε * λ ( p) = −δ λλ one can invert the above to find (3.4) Apart from the consideration of a homogeneous medium, we will also assume that the occupation numbers are independent of color and spin, that is, In order to follow the classical evolution of the system, we need the corresponding equations of motion which read with a current having quadratic and cubic terms in the gauge fields Now we expand the full interacting field according to A a μ = A (0)a μ +δ A a μ , with A (0)a μ a free field and δ A a μ the piece induced by the interactions and given by δ A a μ (x) = − d 4 y iG μν (x − y)J νa (y), (3.8) where we have already used the fact that the propagator is diagonal in color. It is taken to be the retarded one, and in momentum space in the light-cone gauge it reads where → 0 + , while the prescription for the axial pole is irrelevant for our purposes 2 . Now expanding a λa p = a  where we have already dropped the superscript (0) in the field coefficients.
3.1 The Feynman rules for the classical Yang-Mills theory and the diagonal, gain, term Before proceeding to calculate the diagonal contribution to Eq. (3.11) let us establish the Feynman rules for the calculation of δa λa p . Most of the rules remain the same as the corresponding ones in the scalar theory, while we have the modifications listed below.
• Assign a factor V abc μνρ ( p 1 , p 2 , p 3 ) for each cubic vertex and a factor V abcd μνρσ for each quartic one where V abc μνρ ( p 1 , p 2 , p 3 ) = g f abc g μν ( p 1 − p 2 ) ρ +g νρ ( p 2 − p 3 ) μ +g ρμ ( p 3 − p 1 ) ν , (3.12) V abcd μνρσ = −ig 2 f abe f cde (g μρ g νσ − g μσ g νρ ) + f ace f bde (g μν g ρσ − g μσ g νρ ) + f ade f bce (g μν g ρσ − g μρ g νσ ) .  Fig. 7 The contributions to δa p leading to a the loss terms and b the second gain term in the Boltzmann equation in Yang-Mills theory. The gray blob stands for the total amplitude for 2 → 2 scattering stripped off its polarization vectors (denoted by M in the text). In a the momentum k is flowing to the right and the propagator is retarded while in b the momentum is flowing to the left and the propagator is advanced In the above, |M| 2 YM is the scattering amplitude squared, summed over all initial and final colors and polarizations, at order g 4 in the Yang-Mills theory and it reads 4 Fig. 7.b. Putting everything together we find the Boltzmann equation in classical Yang-Mills theory, that is, with |M| 2 YM given earlier in Eq. (3.19). Notice that the combination which appears in the integrand is really the amplitude squared averaged over the color and polarization of the measured gluon and summed over the colors and polarizations of the remaining three gluons. As in the scalar field theory, a factor 1/2 in front of the integral is a symmetry factor due to the fact that particles 2 and 3 are identical.
Before closing let us repeat here an observation made in [21]. At the level of the classical approximation, since f p 1, we can assume a modified definition of the occupation number by replacing f p in the r.h.s. of Eq. (3.5) with f p + 1/2. Then such a replacement is carried over to all occupation numbers appearing in the collision integral of the Boltzmann equation and one sees that the cubic in f terms remain unaltered as they should. Interestingly enough, the generated quadratic in f terms are exactly those present in the more general Boltzmann equation which is valid for all values of f p and is given in Eq. (1.1). However, such a replacement also gives rise to terms linear in f , which are absent from Eq. (1.1).