Shear viscosity of ultrarelativistic Boson systems in the presence of a Bose-Einstein condensate

We calculate for the first time the shear viscosity of ultrarelativistic Boson systems in the presence of a Bose-Einstein condensate (BEC). Two different methods are used. One is the Grad's method of moments and another is the Green-Kubo relation within a kinetic transport approach. In this work we consider a Boson system with isotropic elastic collisions and a gluon system with elastic scatterings described by perturbation QCD (pQCD). The results show that the presence of a BEC lowers the shear viscosity. This effect becomes stronger for the increasing proportion of the BEC in the Boson system and is insensitive to the detail of interactions.


I. INTRODUCTION
The experiments of heavy-ion collisions at the Relativistic Heavy Ion Collider (RHIC) and at the Large Hadron Collider (LHC) [1][2][3] showed strong indications for the formation of a new thermal state of matter composed of quarks and gluons, the quark-gluon plasma (QGP). Theoretically, the color glass condensate (CGC) effective field theory [4] can describe the evolution of the initially freed gluons to a glasma [5][6][7], which is still far from thermal equilibrium. The formation of QGP from this nonequilibrium glasma is a dynamical thermalization process, which is not understood yet from the first principles within QCD. Efforts from recent investigations [8,9] showed that gluons in the glasma can be highly overpopulated. As a fundamental consequence of quantum statistics, a gluon BEC could appear. The gluon condensation, if it occurs, will accelerate the thermalization process and thus, is regarded as a promising mechanism for the fast thermalization of quarks and gluons produced in heavy-ion collisions at RHIC and LHC. Many recent studies have been devoted to the nonequilibrium dynamics and BEC formation within either kinetic approach [10][11][12][13][14][15][16][17][18][19] or classical field theory [20][21][22][23].
Whether or not gluon BEC can be formed in heavy-ion collisions, is still under debate [13,[24][25][26][27]. In this work we will not touch this issue. We assume the existence of a gluon BEC and would like to study its effect on the shear viscosity of gluons.
The shear viscosity is an important quantity manifesting the transport property of the QGP. Extracted from the flow measurements at RHIC and LHC by comparing with the viscous hydrodynamic calculations, one found small numbers of the shear viscosity over the entropy density ratio (η/s) of the QGP [28]. To understand the nature of the small η/s, is one of the motivations for this work.
From kinetic theory we can realize that stronger in- * xuzhe@mail.tsinghua.edu.cn teractions will lead to faster thermalization as well as smaller shear viscosity [29]. When the gluon condensation accelerates thermalization, it is expected that processes corresponding to the gluon condensation will lower the shear viscosity. However, quantitatively it is nontrivial to confirm this expectation, especially for massless Boson systems. In Sec. II we discuss the description of BEC in kinetic theory from binary processes, especially highlight the rate of condensation. The collision rate of processes involving the BEC is infinitely large, while the collision angle is zero. How these extreme processes affect the shear viscosity is one key ingredient of this work. In Sec. III we derive for the first time the shear viscosity of massless bosons in the presence of a BEC analytically by applying the Grad's method of moments [30][31][32][33]. As a first test case, we assume elastic collisions with constant cross sections and isotropic distribution of collision angles. This is customary in kinetic theory calculations. Later in Sec. V we will relax this restriction. As a complementary method, we numerically calculate the shear viscosity in Sec. IV for the same system and with the same interactions by using the Green-Kubo relations within the Boltzmann Approach of MultiParton Scatterings (BAMPS) [34]. The two independent methods agree perfectly and confirm our analytic and numerical calculations. In Sec. V we consider a realistic gluon system, employing binary pQCD cross sections in the presence of a gluon BEC within the numerical framework BAMPS. We will summarize in Sec. VI. The details of the thermodynamic integrals needed in Sec. III are given in Appendix.

II. KINETIC DESCRIPTION OF BOSE-EINSTEIN CONDENSATION
In this section we give a brief description of Bose-Einstein condensation by using the kinetic Boltzmann equation. The detailed description can be found in Ref. [19].
The one-particle phase space distribution function f (x, p) is decomposed into two parts f = f g + f c , where f g denotes the distribution of gas (noncondensate) particles and f c = (2π) 3 n c δ (3) (p) denotes the distribution of the condensate particles with zero momentum. n c (x) is the local particle density of the condensate. In this work we consider elastic collisions only. This assumption is reasonable for systems with a conserved particle number. For a gluon system, however, number-changing processes may have to be taken into account. In Sec. V we will give a short discussion about why we do not consider number-changing processes of gluons in the present work.
Denoting gas particles by g and condensate particles by c, we consider g + g → g + g, g + c → g + g, and g + g → g + c processes. The Boltzmann equations for gas and condensate particles are then given as follows: where dΓ i = d 3 p i /(2E i )/(2π) 3 and f g/c i = f g/c i (x, p i ) with i = 1, 2, 3, 4 is the distribution function of ith particle. Elastic collisions 34 → 12 and 12 → 34 are determined by the collision kernel |M 34→12 | 2 and |M 12→34 | 2 , which are equal.
The rate of the Bose-Einstein condensation can be obtained by the integration of the right-hand side of Eq.
(2) over dΓ 1 in the local rest frame. The details of the integration can be found in Ref. [19]. The final result is The two terms on the right-hand side of Eq. (3) correspond to kinetic processes for the condensation and the evaporation, respectively. E = E 3 + E 4 is the total energy in the collision, P = |p 3 + p 4 | is the total momentum, and s = E 2 − P 2 is the invariant mass. m denotes the particle mass at rest, which is set to zero throughout this paper. We found [19] that in order to describe the condensation of massless bosons with a finite rate, the ratio |M 34→12 | 2 /s at s = 0 should be nonzero and finite. For isotropic collisions, i.e., the distribution of the collision angle is isotropic, we get |M 34→12 | 2 = 32πsσ 22 , where σ 22 is the total cross section. We see that isotropic collisions with finite cross sections can describe the condensation process of massless bosons with a finite rate. At thermal equilibrium, n c does not change with time. The two integrals in the right-hand side of Eq. (3) cancel. More important is that both integrals will be infinitely large, when f g i takes the fully equilibrated Bose-Einstein distribution function, f g i = 1/(e Ei/T − 1), where T is the temperature, because for instance, ∞ m dE 3 f g 3 is logarithmically divergent for m = 0 (which also holds for finite mass). This indicates that both the rates of g + g → g + c and g +c → g +g are infinitely large at equilibrium. Since the main purpose of this work is to calculate the shear viscosity of a Bose gas in the presence of a BEC, one may argue that the infinite collision rates of g + g → g + c and g + c → g + g would naively lead to a vanishing shear viscosity of the gas (noncondensate) particles. However, the collision rate is not the only determinant of the magnitude of the shear viscosity. Another determinant is the collision angle. In processes g + g → g + c and g + c → g + g with massless particles, all the momenta of the noncondensate particles before and after the collision should be parallel. The only change after g + g → g + c and g + c → g + g processes are the magnitudes of the momenta, but not their directions. This corresponds to zero collision angle. One may argue again that those collisions would not contribute to the shear viscosity. We see that the infinite collision rate and zero collision angle are two extremes in the condensation processes. How the counterbalance of these two extremes will affect the shear viscosity is an interesting and nontrivial issue, which will be addressed in the following sections with two different methods. One is analytical from second-order kinetic theory [33], while another is numerical from the Green-Kubo relation within the transport approach BAMPS [34].

III. SHEAR VISCOSITY COEFFICIENT FROM SECOND-ORDER KINETIC THEORY
Relativistic causal dissipative hydrodynamic equations can be derived from the kinetic theory by applying Grad's method of moments [30]. A detailed prescription for the derivation of shear viscosity from the second-order kinetic theory is given in Ref. [33].
In this work we use the general formula of the secondorder shear viscosity derived in Ref. [33] for a special case: a massless Boson system in one-dimensional Bjorken expansion [35]. The local equilibrium distribution function of bosons is the Bose-Einstein distribution function: where u µ (x) is the four-fluid velocity. For a one-dimensional Bjorken expansion, we have u µ = (t, 0, 0, z)/τ with τ = √ t 2 − z 2 [35]. In the local rest frame, u µ = (1, 0, 0, 0). The shear viscosity is a material property. Its value does not depend on the special form of the collective motion of the matter. The assumption of the onedimensional Bjorken expansion will simplify the calculation of the shear viscosity, because in this case the heat flux q µ vanishes in the local rest frame [36]. In addition, for massless systems, the bulk pressure Π becomes zero [36]. Then, the formula of the second-order shear viscosity from Ref. [33] is reduced to π αβ denotes the shear tensor. We express C 0 explicitly as given in Ref. [33]: (6) The integral can be easily performed and we obtain involves interactions between particles by means of the collision term in the Boltzmann equation (1). The viscosity is the measure of the medium response to a small disturbance, which leads the system to deviate slightly from local equilibrium. The one-particle phase-space distribution function f (x, p) has then the following form (the subscripts are omitted), Using the relativistic Grad's 14-moment approximation [30,37] or variational method [31], φ(x.p) is approximated up to the second order of momentum [33,36], Since we have assumed one-dimensional Bjorken expansion, in the local rest frame, the shear tensor takes the form π µν = diag(0, −π/2, −π/2,π). According to Eqs. (7)-(9), P µν will be calculated when putting π µν into φ(x, p). We will see in Appendix that π µν P µν is proportional toπ 2 and thus,π from the numerator and denominator of Eq. (5) cancels out. The shear viscosity does not depend onπ and can be evaluated explicitly, if the matrix elements of particle interactions, which are involved in the collision term C[f i (x, p i )], are given. We now calculate for the first time the second-order shear viscosity of the noncondensate particles in the presence of a Bose-Einstein condensate. We assume isotropic collisions with a constant cross section. A constant cross section means that the total cross section is independent of s. The assumption is useful for testing the numerical framework, and is always interesting, because in kinetic theory constant isotropic cross sections are very often used. With this assumption, some integrals in Eq. (7) can be carried out analytically. The rest has to be computed numerically. We obtain where The result, Eqs. (10) and (11), is new and nontrivial. The details of the integration are presented in Appendix.
The shear viscosity is proportional to the temperature and inversely proportional to the total cross section. For the absence of the condensate (n c = 0) we obtain k = 1.05. This value is smaller than that (k = 1.2) when assuming Boltzmann statistics (neglecting the Bose factors) [38]. This shows that the Bose factors increase the collision rate, which lowers the shear viscosity.
From Eqs. (10) and (11) we realize that the presence of the BEC lowers the shear viscosity. Since the number density of noncondensate particles is n g = ζ[3]T 3 /π 2 , η(n c )/η(n c = 0) only depends on n c /n g . In other words, the larger the proportion of the BEC in the system, the stronger is the effect of the BEC on the shear viscosity. By keeping n g constant (with a fixed temperature) and varying n c , the proportion of the BEC in the system n c /(n c + n g ) will change accordingly. With T = 0.4 GeV and σ 22 = 1 mb we show in Fig. 1 the shear viscosity as a function of n c /(n c + n g ) by the red dashed curve. We see a decreasing shear viscosity when increasing n c . The main result of this work is that we find a finite and nonzero contribution of the BEC (or g + g → g + c and g + c → g + g processes) to the shear viscosity of massless noncondensate particles. Remember [see Eq. (3)] that the rate for the condensation (g + g → g + c) at equilibrium is infinite. The logarithmic divergence comes from the integration Because the same divergence appears in the rate for the evaporation (g+c → g+g), the net rate for the condensation at equilibrium is zero. This shows that the gain and loss term cancel, when the collision term of the Boltzmann equation is integrated in momentum space. In the case of the calculation of the shear viscosity [see Eqs. (5) and (7)], the second moment of the collision term instead of the collision term itself is integrated. The gain and loss terms do not cancel, but lead to a net momentum transfer. In the single gain and loss terms, there are still divergences due to the same reason as in the calculation of the collision rates. We sum all gain and loss terms with divergences and find that the term corresponding to the momentum transfer in processes 1 and there are only integration with the mixed term and no integration such as BE . The divergences in the gain and loss terms cancel with each other. The net momentum transfer is finite and nonzero.
Qualitatively, collinear collisions involving a massless BEC particle (g + g → g + c and g + c → g + g processes) lead to a redistribution of the magnitude of momentum, which is obviously a momentum degradation among different nearby fluid layers, although the collision angle in such collinear collisions is zero. The redistribution of the magnitude of momentum gives a nonzero momentum transfer and a nonzero contribution to the entropy production. Therefore, collinear collisions involving a BEC particle should contribute to the shear viscosity, as demonstrated mathematically in Appendix. In the next section we will use a different method to calculate the shear viscosity in the presence of a BEC, in order to prove the result shown in this section.

IV. SHEAR VISCOSITY COEFFICIENT FROM GREEN-KUBO RELATIONS
According to the several works by Green and Kubo [39,40], which are motivated by Onsager's regression hypothesis [41], transport coefficients can be related to the correlation function of the corresponding flux or tensor in thermal equilibrium. The Green-Kubo relation for the shear viscosity is where i, j = x, y, z and · · · denotes the ensemble average in thermal equilibrium. The sum over i and j gives π ij (r, t)π ij (0, 0) = 10 π xy (r, t)π xy (0, 0) . Different from the one-dimensional Bjorken expansion assumed in the previous section, we consider here a homogeneous static system in global thermal equilibrium. At thermal equilibrium, fluctuations are still present. Thus, π xy (r, t) fluctuates around zero. The dissipation of fluctuations leads to the relaxation of the correlation π xy (r, t)π xy (0, 0) that is determined by the shear viscosity as indicated in Eq. (12).
In this work, fluctuations are realized in BAMPS [42], where test particles are used to represent the particle phase space distribution function. Although the number of test particles should be high enough to ensure the high accuracy of the solution, it is still finite, which then leads to fluctuations of π xy (r, t). A calculation of the shear viscosity using Green-Kubo relations within BAMPS (but without Bose statistics and BEC) has been done in Ref. [34]. We follow this framework, but employ the newly developed BAMPS including Bose statistics and BEC [19]. The most details on the numerical implementations can be found in Refs. [19,34,43]. In the following we briefly present the numerical procedure and show some new numerical implementations.
We consider a static box with volume V = L 3 , where L is the side length. The densities of test particles in the box are the physical particle number densities (n g and n c ) multiplied by N test , which then indicates the number of test particles per real particle. The condensate test particles are approximated as particles with energy being less than ǫ = 2.5 MeV [19]. The chosen cutoff ǫ is small enough, so that this approximation will not destroy the initial equilibrium state [19]. Initially, the spatial coordinates of test particles are sampled homogeneously in the box. While momenta of noncondensate test particles are sampled according to the Bose-Einstein distribution, energies of condensate test particles are sampled according to a uniform distribution within the interval [0, ǫ] and their momenta are isotropically distributed.
The time evolution of test particles in phase space is the consequence of particle collisions and the free moving between two successive collisions. In BAMPS, collisions of particles are simulated in a stochastic way according to collision probabilities corresponding to the collision rates that can be calculated from the collision term in the Boltzmann equation. Collisions are realized in each spatial cell with cell length ∆L within time step ∆t. More details on the numerical implementations of collisions g + g → g + g, g + g → g + c, and g + c → g + g can be found in Ref. [19].
In our calculations we average π xy (r, t) over the space in the box. Thus, at each time t, π xy (t) is evaluated as where the sum is over all N noncondensate test particles in the box. During the time evolution of test particles, π xy will only change, once a collision among test particles occurs, since collisions change the momenta of colliding particles and thus, weakens the correlation π xy (t)π xy (0) . When test particles hit the wall of the box, they will be moved to the opposite wall, so that they can still stay in the box. This periodic boundary condition dose not change π xy , since the momenta of the test particles do not change.
We calculate π xy (t)π xy (0) in one run and obtain the correlation π xy (t)π xy (0) by the average over a large number of runs with different randomly sampled initial conditions. The correlation will be put into Eq. (12) to calculate the shear viscosity. For more technical details refer to Ref. [34].
We present now a new numerical implementation. As we have noticed in Sec. II, the collision rates involving the condensate particle are infinitely large. With the energy cutoff ǫ for the condensate particles, these collision rates are now finite, but still large [19]. In addition, for noncondensate particles, the Bose factor (1 + f g 1 )(1 + f g 2 ) increases the collision rates significantly, when the momenta of colliding particles are small. Figure 2 shows the collision rate per particle with energy E, calculated in one time step. We see that the collision rate increases rapidly, when the energy becomes small. Therefore, on average, particles with smaller energy have larger probability to collide than particles with larger energy within one time step. In BAMPS, we compute the collision probabilities to randomly decide whether a collision occurs.
Since the collision probability is proportional to the time step times the collision rate and the latter becomes huge for particles with small energy, the time step should be chosen quite small, in order to keep the collision probability smaller than 1. This leads to a very time consuming computation, because within one such small time step, particles with large energy have tiny collision probabilities, but have as much numerical operations as particles with small energy. Therefore, the numerical handling for particles with large energy is inefficient. In order to make the computation more efficient, we introduce in this work two kinds of time step, a smaller and a larger one. If the energy of at least one of the two colliding particles is smaller than a cutoff, E cut , the smaller time step is used for calculating the collision probability. If the energies of both colliding particles are larger than E cut , the larger time step is used. In this way we immensely reduce the computing costs. In the calculations we choose E cut = 0.5 GeV empirically.
We list here settings for the numerical calculations. The length of the box is L = 3 fm. The cell length is ∆L = 0.25 fm. We set N test = 2400 and perform 200 independent runs to make ensemble averages. The temperature of the Boson gas is T = 0.4 GeV. We assume isotropic elastic scatterings.
Using the Green-Kubo relation, Eq. (12), the shear viscosity has been calculated for five different cross sections, as shown in Fig. 3. The squares depict the re- sults (multiplied by 10) without a BEC, while the circles depict the results in the presence of a BEC with n c /(n c + n g ) = 30%. The results from the previous section, Eqs. (10) and (11), are shown by the solid and dashed lines for comparisons. We see excellent agreements between the results obtained by using the Green-Kubo relation and the Grad's method of moments.
We have also calculated the shear viscosity for two further BEC proportions, n c /(n c +n g ) = 15%, 50%. In these cases the total cross section is set to be 1 mb. The results together with those for n c /(n c +n g ) = 0%, 30% are shown in Fig. 1 by the circles. Again, we see perfect agreements between the results from two different methods.

V. SHEAR VISCOSITY OF GLUONS IN THE PRESENCE OF A BEC
For a system of massless gluons, the presence of a gluon BEC would lower the shear viscosity as expected from the result in the previous sections. Since the gluon number is not conserved due to number-changing processes such as two gluons go to three gluons and vice versa, whether a gluon BEC can be formed (grow) is still under debates [13,24,25]. Taking a first detailed look at the collision terms corresponding to g + g ↔ g + c + c processes (which are dominant compared to g + g + g ↔ g + c and g + g ↔ g + g + c), we see that is zero for massless gluons in equilibrium, while it is negative for massive gluons with the chemical potential µ = m. This indicates that number-changing processes cannot destroy a massless gluon BEC, while they do for a massive gluon BEC. Following the similar integration that leads to Eq. (3), we find that the Bose-Einstein condensation with a finite rate requires an additional condition: |M 45→123 | 2 /s 2 at s = 4mE − 3m 2 should be nonzero and finite. For m = 0 both the rates of g +g → g +c+c and g +c+c → g +g are infinitely large, similar as the rates of g + g → g + c and g +c → g +g. Therefore, the same numerical handle with the divergence as done for binary collisions in the previous section would be used when including these numberchanging processes. On the other hand, for m > 0, the net rate of the g + g → g + c + c and g + c + c → g + g processes is negative, when µ = m. Its value depends on |M 45→123 | 2 . Assuming |M 45→123 | 2 to be constant, the net rate becomes negative infinite. In this case a BEC may not occur. Whether a pQCD motivated |M 45→123 | 2 will lead to a different result needs more detailed calculations. This and the effect on the shear viscosity will be discussed in a future work. In the current work we assume the dominance of elastic scatterings and ignore number-changing processes g + g ↔ g + c + c, g + g + g ↔ g + c, and g + g ↔ g + g + c. For consistency we also ignore g + g ↔ g + g + g processes, although these processes do not destroy the global equilibrium and could be included to calculate the shear viscosity [34,44,45]. Under these assumptions we expect the (temporarily) presence of a BEC in the early stages of relativistic heavy-ion collisions and calculate in the following the shear viscosity of a gluon gas in the presence of a BEC within BAMPS by using the Green-Kubo relation.
We consider a gluon system in thermal equilibrium with a BEC, which may be formed from a nonequilibrium and overpopulated gluon system produced as an initial state in ultrarelativistic heavy-ion collisions. A simplified form of such an initial distribution of gluons is [8] f init (p) = f 0 θ(Q s − |p|) .
Immediate equilibration would lead to the formation of a BEC for f 0 > 0.154 [10] according to the number and energy conservation. It was shown in Ref. [19] that at equilibrium, where n total = d G f 0 Q 3 s /(6π 2 ) is obtained from the initial distribution (14). d G = 16 is the degeneracy factor of gluons. For Q s = 1 GeV and f 0 = 1.0 we have n c /n total = 37% and T = 0.443 GeV. If gluon scatterings were isotropic, one could obtain the effect of the presence of BEC on the gluon shear viscosity according to Eqs. (10) and (11) : η(37% BEC)/η(No BEC) = 0.846.
Different from isotropic scatterings, the matrix element of the elastic scattering of gluons has the following form: which has been calculated by using the Hard-Thermal-Loop (HTL) treatment [27,46]. s and t are the Mandelstam variables and m D is the Debye screening mass. In thermal equilibrium we have m D = √ 4πα s T [19]. We set α s = 0.3. The total cross section is logarithmically divergent, which is regularized by an upper cutoff of t [19]. Using BAMPS we evolve an equilibrium gluon system with (or without) a BEC in a box. We set T = 0.443 GeV and n c /n total = 37%, (or 0%). The numerical implementations for pQCD scatterings are the same as established in Ref. [19]. We calculate the shear viscosity of gluons from the Green-Kubo relation given in the previous section and obtain η/s = 0.438(with BEC) and η/s = 0.531(No BEC). Here, s denotes the entropy density. At equilibrium, s = d G 2π 2 T 3 /45. The gluon BEC has no contribution to the total entropy. We find that the ratio of the shear viscosity with BEC over the shear viscosity without BEC is 0.825, which is almost the same as that calculated with isotropic scatterings. This suggests that the effect of the presence of BEC on the shear viscosity is insensitive to the detail of scatterings.

VI. CONCLUSIONS
In this paper we have calculated for the first time the shear viscosity of ultrarelativistic Boson systems in the presence of a BEC. For a special case of massless bosons with isotropic elastic collisions, the shear viscosity can be derived analytically by applying the Grad's method, see Eqs. (10) and (11). We found that the presence of BEC or more precisely, the interactions corresponding to the condensation lower the shear viscosity. The larger the proportion of BEC in the Boson system, the stronger is the reduction of the shear viscosity. This analytical result is confirmed by comparing with the numerical result obtained by using the Green-Kubo relations within the transport approach BAMPS. The agreement between these results in turn demonstrated the correct numerical implementations in the BAMPS simulations.
The advantage of the BAMPS simulations over the integrals in the Grad's method is that the computational expenses in the BAMPS simulations are same for any form of the matrix element of elastic scatterings, while the integrals in the Grad's method can be carried out or reduced to integrals with lower dimensions only for simple forms of the matrix element. In addition, the current BAMPS simulations can be easily extended to apply to systems with massive particles. For a gluon system, the matrix element of elastic scatterings calculated from pQCD is more complicated than the isotropic form. We, thus, use the BAMPS simulations to calculate the shear viscosity of a gluon system in the presence of a BEC. Our results showed that a potential formation of a BEC in the early stage of heavy-ion collisions will reduce the shear viscosity of gluons. The reduction of the shear viscosity due to the presence of BEC is insensitive to the detail of interactions. Thus, the analytical formula of the shear viscosity, Eqs. (10) and (11), obtained for isotropic scatterings, can be used to estimate the effect of the BEC on the shear viscosity of gluons.
The sum of integrals with zero power of φ i vanishes, since the collision term vanishes at equilibrium. We calculate integrals with first power of φ i . Integrals with higher power are neglected. Therefore, π µν P µν is proportional toπ 2 .π 2 from π αβ π αβ and π µν P µν in Eq. (5) cancel out. Thus, η does not depend on the magnitude ofπ. Some integrals can be carried out analytically. The rest has to be calculated numerically. The result is where Remember that E = E 3 + E 4 is the total energy and P = p 3 + p 4 is the total momentum. A 1 , A 2 , and B can only be calculated numerically. In the following we carry out explicitly three integrals in Eq. (A1) as examples.

First integral
This integral corresponds to scatterings among noncondensate particles, g + g → g + g. Remember that |M 34→12 | 2 = 32πsσ 22 for isotropic elastic scatterings, where s = (p 1 + p 2 ) 2 = (p 3 + p 4 ) 2 . We have assumed constant cross section σ 22 . At first, we integrate over dΓ 2 with help of δ (3) (p 3 + p 4 − p 1 − p 2 ) and obtain F (p 1 ) indicates the energy conservation Second we integrate over dΓ 1 . To do it, we rotate the coordinate systemp to a new onep ′ withp ′ z paralleling to P. In the new coordinate system, E, |P| ≡ P , E 1 , and the angle between P and p 1 are unchanged. p 1z is transferred to where β x , β y , and β z are the cosine of angles between P and the old coordinate axes,p x ,p y , andp z , respectively. It is easy to prove that β 2 x + β 2 y + β 2 z = 1. The integral over dΓ 1 in Eq. (A8) is evaluated in the new coordinate system where By integrating over φ ′ 1 , Eq. (A11) is equal to ] is the solution of E ′ 1 in F = 0, we perform the integral in Eq. (A13) over E ′ 1 and obtain Putting the integral over Γ 1 , Eq. (A15), and |M 34→12 | 2 = 32πsσ 22 into Eq. (A8), we have where s = E 2 − P 2 , β z = P ·p z /P = (p 3z + p 4z )/P , and Integrating over φ 3 and φ 4 gives where u 3 = cos θ 3 and u 4 = cos θ 4 . By further integral over u 3 , u 4 , E 3 , and E 4 we obtain

Second integral
The only difference of I 2 from I 1 is the additional multiplier f BE 2 (E 2 ), which changed to f BE (E − E 1 ) by the integral over Γ 2 due to the energy conservation. We follow the integration in the previous section until Eq. (A15). Instead of Eq. (A15) we have now which can only be calculated numerically.

Third integral
The only difference of I 3 from I 1 is the additional multiplier f c 2 = (2π) 3 n c δ (3) (p 2 ). Thus, this integral corresponds to scatterings between condensate and noncondensate particles, g + g → c + g or c + g → g + g. Remember that for these scatterings to occur, |M 34→12 | 2 /s should be finite at s = 2mE, see Eq. (3). For isotropic scatterings |M 34→12 | 2 /s = 32πσ 22 is finite for finite cross sections.
We integrate first over dΓ 1 with help of δ (3) (p 3 + p 4 − p 1 − p 2 ) and obtain where m is the rest mass of particles. We will let m to equal zero at the end of the integration. Using the identity the integral over Γ 2 is expressed to The integration over p 2 gives Putting the integral over Γ 2 , Eq. (A27), into Eq. (A23), we have We define G = (E − √ P 2 + m 2 ) 2 − m 2 . According to Eq. (A17) we obtain where A ≡ (E 3 − m)(E 4 − m)/p 3 /p 4 andφ 4,1 andφ 4,2 are two solutions of φ 4 in G = 0, since cosine of α 1 and α 2 = 2π − α 1 are equal. G = 0 also leads to P 2 + m 2 = (E −m) 2 . The integral over φ 4 gives a factor of 2 because of two solutions and the further integral over φ 3 gives a factor of 2π. Then we have I 3 = 1 (2π) 3π 2 C 0 σ 22 n c dp 3 dp 4