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 is not relied on a doubling of the fields, rather it is based on a diagrammatic approach representing the classical solution to the problem.


I. 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] with v p = p/E p the gluon velocity having unit magnitude, F ext a generic external force and C[f ] the collision term or collision integral accounting for the interactions among gluons. Considering only 2 → 2 elastic scattering this collision term reads C[f ] = 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 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 dp ≡ d 3 p (2π) 3 2E p . (1.3) 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 on-shell 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 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 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 transititons [7,17,19,20].
classical analogues 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 is 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  [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 less than 1/α s , the classical field theory simulations should agree with the Boltzmann equation. And 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

II. 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 perturbation 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, 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 it is straightforward to check by performing the integration over k 0 . More precisely, one finds (2.10) 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 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) pp 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 nonvanishing order in λ ∼ g 2 /M 2 . Eventually this translates to imposing 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.
A. 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 be apparent in a while, let us consider the following particular 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 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 both 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) (p + p 1 − p 2 − p 3 ) and after also using momentum conservation in the DA it finally gives a factor δ (3) pp 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 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. 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 diagram 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 say, just for illustrative purposes, that such diagrams arise from the current J(y) expanded to order λ 2 which can be easily 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 diagram 2.a gives 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, a * p 1 a p 2 a * p 4 a p 5 a p 6 a * p → 2δ (3) 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 δ pp , 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) (p + p 1 − p 2 − p 3 ). Then by putting everything together in Eq. (2.15) we arrive at the order λ 2 loss termṡ f p B λ 2 = − 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 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 let 2f Similarly, the diagram 2.b gives 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 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 diagram 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.
FIG. 3. The g 2 contributions to δa p leading to a gain term in the Boltzmann equation.

B. 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 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 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 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 (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 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 diagram 2.b, it is not hard to convince ourselves that we get a contributioṅ f p C g 4 = 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 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 iṡ f p g 4 = 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 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 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ṁ f p A λg 2 = 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 2λ 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 finds 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 = p 1 . Then the propagator bracket in Eq. (2.47) becomes 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 termṡ f p B λg 2 = − 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 2λ 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ṅ f p C λg 2 = 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) 2λ 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)).

III. 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. II 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 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 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 (0)λa p + δa λa p one finds that the change in the field coefficients is given by 3 Finally the occupation numbers evolve in time according to where we have already dropped the superscript (0) in the field coefficients.
A. 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.

B. Loss and gain from the crossed term
In order to derive the loss terms and the second gain term from the crossed contribution in Eq. (3.11) we need to calculate the g 4 contribution to δa λa p . Regarding the loss term, and given the discussion below Eq. (2.52) at the end of Sect. II, it is not hard to understand that the sum of all possible contributing diagrams is given by a λ 2 a 2 p 2 a λ 4 a 4 * p 4 a λ 5 a 5 p 5 a λ 6 a 6 p 6 ε λ * µ (p)ε λ 1 * µ 1 (p 1 )ε λ 2 µ 2 (p 2 )ε λ 4 * µ 4 (p 4 )ε λ 5 µ 5 (p 5 )ε λ 6 µ 6 (p 6 ) M ba 4 a 5 a 6 νµ 4 µ 5 µ 6 (kp 4 ; p 5 p 6 ) G νρ R (k) M ba 2 aa 1 * ρµ 2 µµ 1 (kp 2 ; pp 1 ), (3.21) with k = p 5 + p 6 − p 4 and where M is the total amplitude M given by the sum of Eqs. (3.15) and (3.16), but stripped off its polarization vectors. Because M is imaginary, in the last factor we have let M → − M * and this is the origin of the minus sign in Eq. (3.21). In Fig. 7.a we show the diagrammatic representation of Eq. (3.21).
The steps to be followed now are of course almost identical to those in the scalar theory.
A notable difference is the Lorentz structure of the propagator in Eq. (3.9), which, since the propagator is eventually put on-shell, can be expressed as a sum over polarization vectors. More precisely one can write −g νρ + n ν k ρ + n ρ k ν n · k = ε ν * λ (k) ε ρ λ (k), (3.22) where, as always, a summation over λ is understood. In total we now have eight polarization a λ 5 a 5 p 5 a λ 6 a 6 p 6 a λ a * p → 2δ λ 1 λ 5 δ λ 2 λ 4 δ λ 6 λ δ a 1 a 5 δ a 2 a 4 δ a 6 a δ p 6 p f p f p 1 f p 2 (3.23) and using Eq. (3.18) we arrive at the loss terṁ f p B = − 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) |M| 2 YM 2(N 2 c − 1) It is more than obvious that if we choose the combination a 1 a 2 a * 4 a * 5 a 6 , instead of a * 1 a 2 a * 4 a 5 a 6 , we shall arrive at the second gain term and the respective diagram is shown in Fig. 7.b. Putting everything together we find the Boltzmann equation in classical Yang-Mills theory, that iṡ f p = 1 4E p dp 1 dp 2 dp 3 (2π) 4 δ (4) (∆p) |M| 2 YM 2(N 2 c − 1) 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].