Far-from-equilibrium kinetic dynamics of $\lambda \phi^4$ theory in an expanding universe

We investigate the far-from-equilibrium behavior of the Boltzmann equation for a gas of massless scalar field particles with quartic (tree level) self-interactions ($\lambda \phi^4$) in Friedmann-Lemaitre-Robertson-Walker spacetime. Using a new covariant generating function for the moments of the Boltzmann distribution function, we analytically determine a subset of the spectrum and the corresponding eigenfunctions of the linearized Boltzmann collision operator. We show how the covariant generating function can be also used to find the exact equations for the moments in the full nonlinear regime. Different than the case of a ultrarelativistic gas of hard spheres (where the total cross section is constant), for $\lambda \phi^4$ the fact that the cross section decreases with energy implies that moments of arbitrarily high order directly couple to low order moments. Numerical solutions for the scalar field case are presented and compared to those found for a gas of hard spheres.


I. INTRODUCTION
The Boltzmann equation plays an important role in our understanding of the properties of dilute gases [1]. In the relativistic regime, the Boltzmann equation [2] has been widely applied to describe phenomena in many different fields, ranging from cosmology [3] to heavy-ion collisions [4][5][6][7][8][9][10]. This nonlinear integrodifferential equation describes how the single-particle distribution function evolves in phase-space in the presence of collisions among the constituents of the gas. From the single-particle distribution function quantities associated with conservation laws such as the energy-momentum tensor, or the particle current, can be reconstructed, which provides a way to determine the hydrodynamic evolution of the system.
The properties of the non-relativistic Boltzmann equation have been widely investigated [1]. In fact, besides extensive numerical solutions, an exact solution for a non-relativistic gas of Maxwell molecules has been derived [11][12][13]. This result, known as the Bobylev-Krook-Wu solution, was found by deriving an exact set of coupled nonlinear differential equations for the moments of the distribution function, which then admit an analytical solution for a given choice of initial conditions. Following along the work by Bobylev, Krook, and Wu, in [14,15] the relativistic Boltzmann equation for a gas of particles interacting with constant cross section, in an expanding spacetime, was rewritten in terms of an exact (infinite) set of coupled differential equations for moments of the single-particle distribution function. In that case, because the cross section is independent of the energy, the equations for the moments can be solved recursively, with the solution of moments of order n only depending on moments of order k < n. This property was crucial to find one analytical solution for the moments, which in turn led to the first (and so far, only) analytical solution of the full nonlinear Boltzmann equation for a relativistically expanding gas [14]. This exact solution was then compared to approximate solutions of the Boltzmann equation obtained by employing the relaxation time approximation [16,17] and the linearized collision term, which was useful to determine the validity of such approximations [15].
In this paper we go beyond [14,15] and consider the Boltzmann equation for a gas of classical massless scalar field particles with quartic (tree-level λφ 4 ) self-interactions in a homogeneously expanding, isotropic spacetime. A new covariant generating function method for the moments of the Boltzmann distribution function is introduced in this work to analytically determine the subset of scalar eigenfunctions of the linearized Boltzmann collision operator, and the corresponding eigenvalues. This covariant generating function is then employed to find the exact equations of motions for the scalar moments in the full nonlinear regime. Unlike the case of hard spheres, for λφ 4 the cross section decreases with energy and, as shall be demonstrated in this paper, this implies that moments of all orders are coupled to each order. This coupling prevents finding a full analytical solution, but numerical solutions of these moment equations can be obtained using simple numerical schemes. A comparison between such numerical solutions, and the corresponding solutions found for a gas of hard spheres with the same initial conditions, is presented in this work. This paper is organized as follows. In Section II we discuss the general properties of the Boltzmann equation in Friedmann-Lemaitre-Robertson-Walker spacetime. In Section III we introduce the generating function method and use it to analytically derive the scalar subset of the spectrum and eigenfunctions of the linearized Boltzmann operator for λφ 4 theory. The scalar moment equations of the nonlinear Boltzmann equation for a scalar field are derived in Section IV. Section V includes some numerical results and comparisons between the scalar field and constant cross section solutions. A summary of our findings and the conclusions we have drawn are presented in section VI. Technical details concerning the properties of associated Laguerre polynomials are given in Appendix A. In Appendix B we show how the new generating function method introduced in this paper can be used to derive the exact moment equations for the constant cross section case, which were originally obtained in [14,15] using different methods.
Notation: We use a mostly minus metric and natural units. Four-vectors are defined as a µ = (a 0 , a) and we use · for the scalar product between spatial vectors, i.e., a i b i = a · b.

II. BOLTZMANN EQUATION IN EXPANDING SPACETIME
We consider a homogeneous and isotropically expanding system of massless scalar particles with quartic self-interactions, embedded in a curved spacetime described by the Friedmann-Lemaitre-Robertson-Walker (FLRW) metric [18] (with zero spatial curvature). In the conformal gauge the line element is (1) We note that the FLRW metric g µν written above is related to the standard Minkowski metric via Weyl rescaling [18]. This fact will play an important role when solving the Boltzmann equation in this curved spacetime, as we explain below. Note that in these coordinates the fluid 4-velocity is u µ = (1/a(τ ), 0, 0, 0), and the expanding FLRW geometry induces a nonzero fluid expansion rate θ(τ ) = ∇ µ u µ = ∂ µ ( √ −g u µ )/ √ −g = 3Da/a 2 , where u µ ∇ µ = D and √ −g = a 4 (τ ), with g being the determinant of the metric in (1). Furthermore, for the FLRW metric above there are many nonzero Christoffel symbols, all of them equal to Da(τ )/a(τ ). We take the probe limit in which the energy and momentum of the kinetic particles are negligible in comparison to other sources that define the underlying cosmological scale factor a(τ ) of the metric. In this limit our results are valid for any a(τ ) > 0 (for instance, for a radiation dominated universe, a(τ ) ∼ τ ).
The dynamics of the single-particle distribution function, f (x, k), is given by the relativistic Boltzmann equation in curved space [19][20][21][22][23] For massless particles with momentum k µ , the on-shell condition k µ k µ = 0 implies that k 0 = |k| = k 2 x + k 2 y + k 2 z (we only use covariant momenta). Since the FLRW universe is spatially homogeneous and isotropic, the distribution function must be homogeneous in space and only depend on the momentum via u µ k µ . Thus, we write f (x, k) = f k (τ ) from here on.
The symmetries of the FLRW spacetime strongly constrain the form of the conserved currents of the matter. Due to local momentum isotropy, the viscous shear-stress tensor, energy diffusion, and particle diffusion current vanish exactly. Therefore, for the massless gas one may write the energy-momentum tensor, T µν , as and the particle 4-current, N µ , as Above, we introduced the spatial projector orthogonal to the 4-velocity, ∆ µν ≡ g µν − u µ u ν . In FLRW, the total energy ε = T µν u µ u ν and particle n = u µ N µ densities in the local rest frame are only functions of τ . The time evolution of these quantities is fully determined by the conservation laws With initial condition a(τ 0 ) = 1, they are solved by n(τ ) = n 0 /a 3 (τ ) and ε(τ ) = ε 0 /a 4 (τ ), where n 0 and ε 0 are constants. Furthermore, for a conformal gas, one may write ε ∼ T 4 and n ∼ T 3 in terms of a suitably defined temperature T , such that T (τ ) = T 0 /a(τ ) and T 0 is the initial temperature scale.
As discussed in detail in [22], for a conformal system the Boltzmann equation transforms covariantly under a Weyl transformation of the metric g µν → e −2Ω g µν . This is the case for the massless (on-shell) λφ 4 theory considered in this paper, where the corresponding interaction cross section does not break conformal invariance. This property implies that, if one conveniently writes the metric in conformal gauge, the factors of a(τ ) will cancel and one can just solve the Boltzmann equation on the flat piece of (1). We will use this property to perform all of our calculations in Minkowski spacetime, leading to a fluid 4-velocity that is static, u µ = (1, 0, 0, 0). To recover the time dependence effect from the metric one can just use the well-known rules based on Weyl rescaling to obtain the quantities at hand (e.g., if ε 0 is the energy density computed using the flat dynamics, the energy density after recovering the Weyl factor will be simply ε 0 /a(τ ) 4 ) [21,22]. Another useful consequence is that the left-hand side of (2) considerably simplifies and the Boltzmann equation becomes where we defined the scalar E k = u µ k µ . Also, after the Weyl rescaling of the metric, we rescale all the momenta (k µ → k µ /T 0 ) and time (τ → τ T 0 ) so both sides of the Boltzmann equation written above are dimensionless. We will work with those rescaled (dimensionless) quantities throughout the paper, unless otherwise specified. In this work we will only consider the classical limit (Boltzmann statistics) where the Boltzmann equation for on-shell massless scalar particles with a λφ 4 interaction (i.e., a cross section σ(s) ∼ λ/s, where s is the square of the center of mass energy) can be written as [2] where λ is a dimensionless constant that denotes the strength of the collisions, and is the Lorentz invariant momentum space integral [2]. We note that equilibrium distribution function is given by where α > 0 is the fugacity, and this function is a zero of the collision term, i.e. C[f eq k ] = 0 [2]. In the next section we develop a covariant generating function method and use it to determine the exact set of eigenvalues and eigenfunctions of the scalar part of the spectrum of the linearized collision operator associated with (7). We then show in Section IV that the method can also be used to determine the exact set of equations of motion for suitably defined scalar moments of the distribution function, which we solve to determine the corresponding solution of the Boltzmann equation.

III. SCALAR SPECTRUM OF THE LINEARIZED COLLISION OPERATOR
The linearized Boltzmann equation is framed in terms of perturbations about the equilibrium distribution function, truncated at first order in deviations. The distribution function is written as where φ k parameterizes the deviations from equilibrium. After substituting this expression into the Boltzmann equation (7), the result to linear order in φ k is Defining the linearized Boltzmann collision operator as the linearized Boltzmann equation can then be written as where we have decomposed the linearized collision operator into a gain term and a loss term which are respectively given by To determine the spectrum of the linearized Boltzmann operator, we propose an Ansatz for the eigenfunctions in terms of the associated Laguerre polynomial, L (1) n (E k ) [24]. We note that as shown in [14], so The first integral can be evaluated explicitly to obtain k f eq k = α 2π 2 . The second integral must be zero because of the orthogonality of the associated Laguerre polynomials. Therefore, the loss term becomes This result indicates that, if L (1) n is an eigenfunction of L, it will be an eigenfunction of both L loss and L gain , separately.
The gain term is considerably more difficult to investigate. For this purpose, it is convenient to use the generating function for the associated Laguerre polynomials [see (A3)] with v ∈ [0, 1). Further properties of this generating function and the associated Laguerre polynomials are discussed in Appendix A. Motivated by seminal work of Refs. [11][12][13], we define the following quantity which can also be written as where we defined the total 4-momentum P µ T = k µ + k µ . Since the integral over pp is a Lorentz scalar that depends only on the timelike 4-vectors u µ and P µ T , the result of these integrals can only depend on those 4-vectors via u µ P µ T and s = P T µ P µ T , with s being the traditional Mandelstam variable [2]. This implies that this integral is invariant under the exchange of u µ and the normalized total 4-momentum,P µ T ≡ P µ T / √ s. This interchange leads to The integral I is a Lorentz scalar and, therefore, can be calculated in any Lorentz frame. For the sake of convenience, we perform this task in the local rest frame of the system where u µ = (1, 0, 0, 0). In this frame, Then, the integral simplifies to wherep = p/E p = p/p is a unit vector in three dimensions. We further define cos θ kp = x and cos θ k p = y, where θ kp is the angle between k and p, and θ k p is the angle between k and p, respectively. This allows I to be written as This integral can now be evaluated using standard techniques to find that Expanding term by term in powers of v, the gain term is then given by The gain term and the loss term are then combined to obtain where the eigenvalues are given by Therefore, we see that L (1) n (E k ) and χ n are, respectively, the exact eigenfunctions and eigenvalues of the scalar part of the spectrum of the collision operator for a massless gas of particles with quartic self-interactions. As n → ∞, the eigenvalues approach − αλ 4π 2 , which indicates that the lifetime of the higher-order modes approaches 4π 2 αλ . Finally, we note that going back to standard units where the momenta are not scaled by T 0 , one finds L (1) n (E k /T 0 ) and that the eigenvalue χ n has dimensions of T 2 0 .

IV. EXACT EQUATIONS OF MOTION FOR THE MOMENTS
The goal of this section is to rewrite the full Boltzmann equation in (7) in terms of ordinary differential equations for suitably defined Lorentz scalar moments of f k , which can be solved using standard numerical routines. Having in mind the results from the previous section, it is natural to expand the distribution function in terms of an associated Laguerre basis L (β) n (E k ) are the eigenfunctions of the linearized collision operator, it turns out that in order to find the exact equations of motion for the nonlinear case, it is best to consider a basis in terms of L (2) n (E k ). This can be understood as follows. Consider the Boltzmann equation (7) and multiply it on both sides by L (β) n (E k ) (with arbitrary β) and then integrate it over k, which gives While the right-hand side can be simplified using the techniques discussed in the previous section for any integer value of β, a bad choice for the coefficient β above can make the left-hand side unnecessarily complex. Indeed, if β = 1, after performing the angular integrals the left-hand side becomes (apart from constant multiplicative factors) One could then decompose the distribution function in terms of generic associated Laguerre polynomials as follows which can be used back in (31) to find that the left-hand side becomes To avoid having a sum of terms already on the left-hand side of the equations for the moments, it is clear that one should set β = γ = 2. Any other choice would severely complicate our analysis. In fact, in this case the orthogonality condition for the associated Laguerre polynomials (A2) can be used and the moment equations become simply where we defined c n ≡ c (2) n to ease the notation. In this way, the moments are finally defined as in [15] c n (τ ) = 2 (n + 1)(n + 2)n 0 k and, once (34) is solved and the c n moments are found, one can always recover back the distribution function as follows: We note that, due to the conservation laws, c 0 = 1 and c 1 = 0 at all times. We will now proceed to express the right-hand side of (34) directly in terms of the c n moments. The result of the integral can be found using the generating function of the associated Laguerre polynomials as follows. We define where we introduced generating function related to the gain and loss terms of the collision term, in such a way that This provides a way to determine the integral in (37). The loss term is simpler and, thus, is evaluated first. Using the trivial identity the loss term can be re-written as These exponentials can be written in terms of associated Laguerre polynomials using (for a ≥ 0) so the loss term becomes We then evaluate the integrals over a, b and use (35) to write the loss term as To compute the gain term, we first evaluate the integral where we have once again used the symmetries of the integral to switch u µ ↔P µ T , as was done when deriving the eigenfunctions of the linearized collision operator. We thus arrive at The gain term is then given by As in the linearized case, we define x = cos θ kp and y = cos θ k p so that the gain term becomes By defining X = 1 , and using (42), we obtain Once again using (44), this can be expressed as Evaluating the integrals over a, b, X, Y using standard techniques and combining this with (35), we find that the gain term is given by All that remains is to combine the gain and loss terms. Once that is done, one finds After performing several redefinitions of summation indices and using the expansion valid for |v| < 1, we obtain This has the same form as (41) except it is missing the first two terms, N = 0, 1, which are fixed by the conservation laws. So, we use the generating function to determine that where the first two terms will not contribute because Dc 0 = Dc 1 = 0. Comparing this to (56) we arrive at a set of equations for the evolution of the moments, for N ≥ 2. This defines a set of equations that determines the exact time evolution of the moments c N , with which a full solution for the distribution function of the Boltzmann equation can be reconstructed using Eq. (36). As such, the equations above can be used to determine how an arbitrarily far from equilibrium state of the gas of massless scalar particles evolves in time in an expanding FLRW universe. To the best of our knowledge, this is the first time the exact set of equations of motion for the (scalar) moments of the full nonlinear Boltzmann equation describing massless scalar particles has been derived. It is instructive to compare our result in (58) to the corresponding set of equations derived in [14,15] that describes a massless gas with constant cross section in FLRW. In the latter, because of the drastic assumption about the particle interactions, the solution of the n-th moment only depended on the dynamics of the previous moments. Therefore, for the system considered in [14,15] an iterative procedure could be easily employed to obtain the moments for arbitrary initial conditions. Furthermore, the equations simplified so much in that case that even an analytical solution for the moments (and, consequently, to f k ) could be found [14]. In contrast, in the case of λφ 4 interactions where the cross section varies with the energy, the results of this section show that one is still able to find the exact set of equations (58) that describes the evolution of the moments, but now we see that the derivative of the n-th moment does not depend only on the previous moments. Rather, it depends on the sum over all the moments via ∞ l=0 c l . This means that analytical solutions will be even harder to find than before. Furthermore, (58) cannot be solved using a simple iterative scheme. However, (58) can still be solved numerically, as we show in Section V.

Uniqueness of equilibrium
In this section we show that the asymptotic equilibrium state solution of (58) is unique, as expected. The argument works as follows. First assume that all the moments have reached their asymptotic state such that dc N /dτ = 0 for all N . Then prove using the equations of motion (58) that this implies that all moments c N for N ≥ 2 are equal to zero -which indicates that the standard equilibrium state, described by the distribution f eq k , has been reached.
for all N . Above, we introduced the quantity The binomial in the summation above is zero whenever 2 > N − n, so we start by considering the case where N = 2. Then, the only term that will contribute to the sum is that of n = 0, which gives the condition The next case to consider is that of N = 3. Now, two terms will contribute to the sum, those of n = 0, 1. This will then give the condition Using that the values of the two first moments are already known, c 0 = 1 and c 1 = 0, this condition provides an initialization for a proof by induction.
Proceeding with standard induction, it is assumed that c 2 , ..., c k = 0 for some k ≥ 2. Then, N = 1 + k is considered. In this case, we have Since we have assumed that the only non-zero c n for n < k is c 0 = 1, it follows that the term c n−m c m is only non-zero when n − m, m are some combination of 0, k + 1. We also know that the binomial is zero whenever 1 > k − n + 1, which is equivalent to n > k. So, it follows that n can never reach k + 1 and the only non-zero terms will occur when n = m = 0. But, as before this will always cancel with the case of c n in the next term. So, all that remains in the summation is However, the binomial is equal to zero for n > k, so we can really truncate this series at This leaves only one term, By induction, it follows that all c n are equal to zero except for c 0 = 1. This means that the only steady-state solution is the standard equilibrium state, as expected. This provides an alternative way to show that the equilibrium state is unique, starting from the exact moment equations.

A comment on the thermalization time
In the following, we will show that M only vanishes in equilibrium, i.e., when all moments c m≥2 = 0. First, we assume that M is equal to zero for all τ ≥ τ 0 . Using this assumption, the equation of motion for the second moment obtained from (58), evaluated at τ ≥ τ 0 , implies that dc 2 /dτ = 0 and, consequently, that c 2 is a constant. If we then analyze the equation of motion for the third moment, we see that dc 3 /dτ only depends on c 2 and, therefore, this derivative is a constant. This would imply that c 3 is a linear function of time. Similarly, if one analyses the equation of motion for the fourth moment, it will depend only on c 2 and c 3 and, thus, c 4 will be a quadratic function of time. The same argument can be applied to the remaining moments and one will conclude that the moments will all display a polynomial behaviour with time, with the highest possible power of time being τ n−2 for c n . Therefore, this type of solution diverges as τ → ∞, which is incompatible with the existence of the equilibrium state and the fact that moments of f k must be finite. The only solution that does not display this divergent behavior corresponds to the case where c n≥2 = 0the equilibrium solution. Therefore, one can see that the thermalization time of the system is well estimated by the time it takes for M to vanish.

V. NUMERICAL RESULTS AND COMPARISON
Since an exact analytical solution has not been obtained, simple numerical procedures are required to solve the moment equations and compare their solutions to those found in the constant cross section case previously considered in [14,15]. The evolution equation for the moments in the scalar field case contains an infinite summation, which makes it harder to solve in comparison to the constant cross section case. Here solutions are obtained by considering some maximum moment defined by N max and numerically solving the evolution equations up to that highest moment. To choose this maximum moment, it is examined at which point the evolution of the moments and the distribution function reach a steady state and do not appreciably change anymore. For each of the figures presented in this work, the maximum moment N max = 90 is used, and the time is scaled by λn 0 . To solve the differential equation after this cutoff is applied, a fourth-order Runge-Kutta algorithm is applied.
Our results are compared to the constant cross section case so a brief summary of the results found in the latter is included here. In [14,15] an exact differential equation for the moments was derived for the case of massless particles interacting with a constant cross section. However, the differential equation for the n-th moment only depends on moments with order less than or equal to n. This allows for a simple iterative approach to be employed to obtain an exact analytical solution for each of the moments and, therefore, for the full distribution function. In particular, for the initial condition it is found that the corresponding (Laguerre-based) moments for hard spheres evolve as [15] c n (τ ) = 1 − n 4 n e −nτ /6 .
This solution is then used as the basis for numerical simulations of the scalar field solution. In fact, we use (67) as the initial conditions for the moments satisfying the equations of motion (58), derived for the scalar field case. The evolution of the moments and the corresponding distribution function for the scalar field system are then compared to the analytical solution in (68). Figure 1 shows the evolution of the Laguerre moments as a function of time for the scalar field and constant cross section cases. In the constant cross section case, each moment directly approaches the equilibrium value at an exponential rate according to the analytical solution. However, in the scalar field case there is some oscillation in many of the moments. In Fig. 2, the effects of these differences on the distribution function can be seen. In the constant cross section case, the system equilibrates considerably faster, and also more uniformly in time. On the other hand, in the scalar field case, the system initially equilibrates faster for low values of k, before smoothing out as the system approaches equilibrium over time. Figure 3 shows how the scalar field distribution function evolves for early times. It can be seen clearly that the initial equilibration process is much more rapid near zero momentum. It would be interesting to see how this is modified when quantum statistics effects are taken into account.
In Fig. 4 the evolution of 1 + M = ∞ n=0 c n is examined for both interactions. We remind the reader that it was argued earlier in Section IV 2 that this quantity provides an estimate of the thermalization time (and rate). In the constant cross section case, it appears that the rate at which the system approaches equilibrium increases rapidly after early times, until the system almost reaches equilibrium. On the other hand, the rate at which the scalar field system approaches equilibrium initially is much slower, steadily increasing over time before slowing down again as equilibrium is approached. This discrepancy may be explained by the oscillatory behavior of moments seen in Fig. 5. Initially some of the moments actually rapidly diverge from the equilibrium value, before more slowly oscillating back towards the equilibrium value. In other words, the scalar field moments do not monotonically approach equilibrium, unlike the constant cross section where all of the moments exponentially decay to equilibrium.

VI. CONCLUSIONS
In this paper we investigated the dynamics of a gas of massless scalar particles with quartic (tree level) self-interactions in Friedmann-Lemaitre-Robertson-Walker spacetime, described by the Boltzmann equation. We demonstrated that the nontrivial far-from-equilibrium dynamics of this system can be determined by solving an infinite set of ordinary differential equations for suitably defined moments of the distribution function. Unlike the case of a gas with constant cross section considered in Refs. [14,15], the fact that in λφ 4 the cross section σ(s) ∼ 1/s makes deriving the exact set of equations of motion for the moments a much more complex task.
We have overcome this challenge by using a new covariant generating function for the scalar moments of the Boltzmann distribution function (see III and IV). This method, which is a relativistic generalization of the generating function techniques used in [11][12][13], was used here to analytically determine for the first time the eigenfunctions and eigenvalues of the scalar part of the spectrum of the linearized collision operator. The spectrum of the collision operator has never been determined in the relativistic regime and even results in a given subspace (such as the scalar sector) are extremely rare -the only other known result can be found in Ref. [14,15], also in the scalar sector. Furthermore, this covariant generating function was also employed to find, for the first time, the exact nonlinear set of equations of motion for the scalar moments in the full nonlinear regime. We showed that the dependence of the cross section with the center of mass energy implies that moments of arbitrarily high order directly couple to low order moments. This should be compared to the constant cross case studied in [14,15] where the n−order moment only coupled to moments of order m < n.
Numerical solutions for the scalar field case were presented and compared to those found for a gas of hard spheres, for the same set of far-from-equilibrium initial conditions. Overall, we found that the dependence of the cross section with the energy introduces more structure in the time evolution of the system, with the moments in the scalar field case displaying oscillations in their approach to equilibrium, while for the constant cross section example the same moments just quickly exponentially decay towards their equilibrium values. This different behavior has consequences to the distribution function as well, which is reconstructed using the moments. In the gas with constant cross section, the system equilibrates much more quickly and also more uniformly in time. On the other hand, for the scalar field case, the system initially equilibrates much faster for low values of k, before smoothing out as it approaches equilibrium over time. We remark that we only considered classical (Boltzmann) statistics in this work. It would be very interesting (and challenging) to generalize our approach to the case where the bosonic nature of the particles is taken into account, as done in Ref. [25]. In that context, one could investigate if our approach would be useful in the investigation of the far-from-equilibrium dynamics of a system that can Bose condense [26][27][28][29].
Concerning the results of this paper, a clear next step would be to see if our generating function method could be used to derive not only the scalar part of the spectrum, but rather the full set of eigenvalues and eigenfunctions of the linearized collision operator considered here. The breaking of isotropy makes this task much more complex, given that now the generating function would have to produce all the scalar, vector, and tensor sectors of the spectrum.
Using this result, the gain term can now be expressed as Above, x and y have been defined as cos θ kp = x and cos θ k p = y, just as in the case of a scalar field. From here, the integral can be simplified using the change of variables X = 1 2 v 1−v (1 + x) and Y = 1 2 v 1−v (1 + y) and by substituting an expansion of the exponential in terms of the associated Laguerre polynomials from (44). The integrals can then be evaluated either using standard techniques to express the gain term as a sum over the moments Next, we consider the loss term, which is given by where the same definition for x and y are used as in the gain term. These integrals can be evaluated explicitly, except for the integral over k. This gives This integral is evaluated by expanding the exponential as in (44), which gives J loss = λn 2 0 2 ∞ n=0 v n (n + 1)(n + 2)c n .
The gain and loss terms can then be combined and compared term-by-term to obtain a differential equation for the c n moments (defined in (35) where the time τ has been scaled by the mean free path l 0 = 1/n 0 λ. Note that this equation implies that the derivative of the n-th moment depends only on the moments c i for i < n, so the system of equations can be easily solved iteratively.
As was shown in [14,15], for the initial condition In this case each of the moments (with n > 1) directly approaches the equilibrium value of zero rather than exhibiting some oscillatory behavior, as found in the scalar field case. This solution is the simplest avenue to compare the constant cross section results to the corresponding results for the scalar field case using the same initial conditions (67), as done in the main text.