Angular mode expansion of the Boltzmann equation in the small-angle approximation

TL;DR

This paper introduces an angular mode expansion method to solve the Boltzmann equation for a gluon plasma, demonstrating rapid convergence and improved computational efficiency, and investigates thermalization dynamics during longitudinal expansion.

Contribution

The paper presents a novel angular mode expansion approach for the Boltzmann equation, showing its effectiveness and faster convergence compared to traditional methods in expanding gluon plasma scenarios.

Findings

Mode expansion converges rapidly for practical cases.

Thermalization involves pressure equilibration lag due to expansion.

First angular mode relaxation is hindered by expansion effects.

Abstract

We use an expansion in angular mode functions in order to solve the Boltzmann equation for a gluon plasma undergoing longitudinal expansion. By comparing with the exact solution obtained numerically by other means we show that the expansion in mode functions converges rapidly for all cases of practical interest, and represents a substantial gain in numerical effort as compared to more standard methods. We contrast the cases of a non expanding plasma and of longitudinally expanding plasmas, and follow in both cases the evolutions towards thermalization. In the latter case, we observe that, although the spherical mode function appears to be well reproduced after some time by a local equilibrium distribution function depending on slowly varying temperature and chemical potential, thereby suggesting thermalization of the system, the longitudinal and transverse pressures take more time to equilibrate. This is because the expansion hinders the relaxation of the first angular mode function. This feature was also observed in a simpler context where the Boltzmann equation is solved in terms of special moments within the relaxation time approximation, and attributed there to the particular coupling between the first two moments of the distribution function. The present analysis confirms this observation in a more realistic setting.

Equations

$$\left(\frac{\partial }{\partial \tau }-\frac{p_z}{\tau }\frac{\partial }{\partial p_z}\right)f(\tau ,\mathbf{p})=-C[f],$$

$$p\equiv |\mathbf{p}|\text{and}\cos \theta \equiv \frac{p_z}{p}.$$

$${\mathcal{L}}_n(\tau )=\int \frac{d^{3}p}{(2\pi {)}^3}pf(\tau ,\mathbf{p})P_{2n}(\cos \theta ),$$

$$f(\tau ,p,\cos \theta )=\sum _{n=0}^{\infty }{f}_n(\tau ,p)P_{2n}(\cos \theta ).$$

$$f_n(\tau ,p)=\frac{4n+1}{2}{\int }_{-1}^{1}dxf(\tau ,p,x)P_{2n}(x).$$

$${\mathcal{L}}_n(\tau )=\frac{2}{4n+1}\frac{1}{4{\pi }^2}{\int }_0^{\infty }dp{p}^3{f}_n(\tau ,p).$$

$$\mathcal{N}=\frac{{\nu }_{\text{g}}}{2{\pi }^2}{\int }_0^{\infty }dp{p}^2{f}_0(\tau ,p),$$

$$\mathcal{E}=\frac{{\nu }_{\text{g}}}{2{\pi }^2}{\int }_0^{\infty }dp{p}^3{f}_0(\tau ,p),$$

$${\mathcal{P}}_L={\nu }_{\text{g}}\frac{{\mathcal{L}}_0+2{\mathcal{L}}_1}{3},{\mathcal{P}}_T={\nu }_{\text{g}}\frac{{\mathcal{L}}_0-{\mathcal{L}}_1}{3}.$$

$$C[f(\tau ,\mathbf{p})]=-4\pi N_c^2{\alpha }_s^2{l}_{\mathrm{Cb}}{\nabla }_{\mathbf{p}}\cdot \left[I_a{\nabla }_{\mathbf{p}}f(\tau ,\mathbf{p})+I_b\frac{\mathbf{p}}{p}f(\tau ,\mathbf{p})(1+f(\tau ,\mathbf{p}))\right],$$

$$I_a(\tau )=\int \frac{d^{3}p}{(2\pi {)}^3}f(\tau ,\mathbf{p})\left[1+f(\tau ,\mathbf{p})\right]$$

$$I_b(\tau )=2\int \frac{d^{3}p}{(2\pi {)}^3}\frac{1}{p}f(\tau ,\mathbf{p}).$$

$$I_b(\tau )=\frac{1}{{\pi }^2}{\int }_0^{\infty }dpp{f}_0(\tau ,p).$$

$$I_a(\tau )=\frac{1}{2{\pi }^2}{\int }_0^{\infty }dp{p}^2\left\{f_0(\tau ,p)+\sum _{n=0}^{\infty }\frac{1}{4n+1}{\left[f_n(\tau ,p)\right]}^2\right\}.$$

$$\frac{\partial }{\partial \tilde{\tau }}{f}_n(p)$$

$$-C_n[f],$$

$$C_n[f(p)]$$

$$-I_b\frac{1}{p^2}\frac{\partial }{\partial p}{p}^2\left[f_n(p)+\frac{4n+1}{2}\sum _{m,l=0}^{\infty }{A}_{n,m,l}{f}_m(p)f_l(p)\right],$$

$$\tilde{\tau }=4\pi N_c^2{\alpha }_s^2{l}_{\mathrm{Cb}}\tau .$$

$$a_n=\frac{8n^2+4n-1}{(4n-1)(4n+3)},b_n=\frac{(2n-1)2n}{(4n-3)(4n-1)},c_n=\frac{(2n+1)(2n+2)}{(4n+3)(4n+5)},$$

$${\tilde{a}}_n=\frac{2n(2n+1)}{(4n-1)(4n+3)},{\tilde{b}}_n=-\frac{(2n-2)(2n-1)2n}{(4n-3)(4n-1)},{\tilde{c}}_n=\frac{(2n+1)(2n+2)(2n+3)}{(4n+3)(4n+5)}.$$

$$A_{n,m,l}={\int }_{-1}^{1}dx{P}_{2n}(x)P_{2m}(x)P_{2l}(x).$$

$$\frac{\partial \mathcal{N}}{\partial \tau }=-\frac{\mathcal{N}}{\tau }.$$

$$\frac{\partial \mathcal{E}}{\partial \tau }=-\frac{\mathcal{E}+{\mathcal{P}}_L}{\tau }$$

$$\frac{\partial \mathcal{E}}{\partial \tau }=-\frac{4}{3}\frac{\mathcal{E}}{\tau }(\text{for }N=0).$$

$$f({\tau }_0,\mathbf{p})$$

$$f_n({\tau }_0,p)$$

$$=f_{\mathrm{ini}}{\delta }_{n,0}\Theta (Q_s/{\xi }_0-p)$$

$$+f_{\mathrm{ini}}\left[P_{2n+1}\left(\sqrt{\frac{Q_s^2-p^2}{({\xi }_0^2-1)p^2}}\right)-(1-{\delta }_{n,0})P_{2n-1}\left(\sqrt{\frac{Q_s^2-p^2}{({\xi }_0^2-1)p^2}}\right)\right]\Theta (Q_s-p)\Theta (p-Q_s/{\xi }_0).$$

$$\frac{1}{{\tau }_R}\sim 4\pi N_c^2{\alpha }_s^2{l}_{\mathrm{Cb}}\frac{I_a}{p^2}.$$

$$I_a(\tau )\simeq \frac{\mathcal{N}}{{\nu }_{\text{g}}}=\frac{Q_s^3{f}_{\mathrm{ini}}}{6{\pi }^2{\xi }_0}\frac{1}{Q_s\tau }.$$

$${\tau }_R(\tau )\sim \frac{3\pi {\xi }_0}{2N_c^2{\alpha }_s^2{l}_{\mathrm{Cb}}{f}_{\mathrm{ini}}}\tau .$$

Full text

• a.

Institut de Physique Théorique, Université Paris-Saclay

CEA, CNRS, F-91191 Gif-sur-Yvette, France

• b.

European Centre for Theoretical Studies in Nuclear Physics

and Related Areas (ECT* ) and Fondazione Bruno Kessler

Villa Tambosi, Strada delle Tabarelle 286, I-38123 Villazzano, Italy

1 Introduction

This paper concerns the kinetic description of expanding quark-gluon plasmas, an area where there has been a lot of activity over the last years. The goals of these calculations is to understand the behavior of matter produced in ultra-relativistic heavy ion collisions, in particular the emergence of hydrodynamics at late time. As representative of some recent works on the subject, let us mention those which concern the solution of the Boltzmann equation with elastic and inelastic processes (see e.g. [1, 2, 3]), or the recent work focusing on very small coupling, of very high occupancies where simple scaling laws emerge [4]. Aside from such elaborate simulations which aim at a realistic description of physical systems, while incorporating detailed QCD contributions to the collision kernel, there is a need to simplify the theoretical description, and get under good theoretical control the various important physical phenomena that are involved, and their associated time scales. The present work is an effort in this direction.

We shall consider a simple setting, that of an expanding plasma composed of gluons only, and with Bjorken symmetry [5]. Our main concern will be to understand how the momentum distribution gets isotropic as a result of the collisions, and this in spite of the strong longitudinal expansion. One of our goal here is to extend the analysis carried out in [6, 7], where it has been pointed out that, within kinetic theory, the competition between free streaming and collisions could be understood in terms of equations for simple moments of the distribution function. These moments, which we called ${\cal L}$-moments, were introduced in [8]. Their knowledge is not sufficient to fully reconstruct the distribution functions, but they capture the essential features of the angular dynamics. In the case where the collision kernel of the kinetic equation is given by a relaxation time approximation, the equations for the ${\cal L}$-moments form an infinite hierarchy that can be solved by truncation [6]. Such a procedure appears to converge quickly, and even the lowest order truncation that takes into account only the first non trivial distortion of the momentum distribution suffices to get an accurate description. We believe that this particular structure is robust, and part of our purpose in this paper is to show that this is indeed the case, by examining a more realistic collision kernel, that of elastic QCD interactions within the small scattering angle approximation [9, 10, 11].

When such a collision kernel is taken into account, the equations for the ${\cal L}$-moments do not close anymore, as new types of moments appear. This is because the ${\cal L}$-moments do not fully take into account the radial momentum dependence (only the angular part is treated exactly, the radial momentum being treated in an average way). There exists of course a number of methods based on more complete moments that could be envisaged (see e.g. [12, 13]). However, we want to retain the simplicity of the angular structure that emerged in our previous analysis, and explore another approach in which the distribution function is expanded in angular mode functions $f_{n}(p)$. The basic assumption is that only a few terms in this expansion are needed, as was the case in the expansion in ${\cal L}$-moments. We shall see that this is indeed the case, which brings a significant simplification in the solution of the Boltzmann equation, since, for the case considered, it reduces a two-dimensional equation to a finite (small) set of one-dimensional equations.

The outline of this paper is as follows. In the next section, we describe the angular mode expansion that we shall be using, and present the resulting form of the kinetic equation which turns into an infinite hierarchy of coupled equations. In the following section, we discuss the initial conditions that are used in the numerical calculations, and we give an estimate of the relaxation time scale, or equivalently of the collision rate. The qualitative features of the evolution of the system in the expanding case depend on the ratio of this collision rate to the expansion rate. We then present the results of a set of numerical calculations. First we consider the non expanding case, where we can analyze in details the isotropization of the momentum distribution and the approach to equilibrium. Then we turn to the case of a longitudinally expanding system, where the approach to equilibrium is studied as a function of the ratio between the collision rate and the expansion rate. Some technical material is presented in two Appendices.

2 Angular mode expansion

The Boltzmann equation for a longitudinally expanding system takes the form

[TABLE]

where $f(\tau,{\boldsymbol{p}})$ denotes the single particle distribution function and $C[f]$ is the collision kernel. We have assumed here that the system is boost-invariant in the longitudinal $z$-direction, and uniform in the transverse $(x,y)$-plane. The distribution function $f$ depends then only on the proper time $\tau=\sqrt{t^{2}-z^{2}}$ and on the momentum ${\boldsymbol{p}}=(p_{x},p_{y},p_{z})$. In writing Eq. (1), we have also exploited boost invariance, and Eq. (1) holds in the central slice around $z=0$.

We further assume that the momentum distribution function is isotropic in the ${\boldsymbol{p}}_{\perp}=(p_{x},p_{y})$-plane. In this case, we can express it as a function of

[TABLE]

In Eq. (1), the drift term $-\frac{p_{z}}{\tau}\frac{\partial}{\partial p_{z}}f$ is responsible for the effects of the longitudinal expansion which tends to make the momentum distribution anisotropic, while the collision term tends to make it isotropic. Our primary interest in this paper is to study the competition of these two effects, and follow the evolution towards isotropization and hydrodynamic behavior. To do so, we generalize the approach initiated in Refs. [8, 6] which is based on the following moments of the distribution function

[TABLE]

where $P_{2n}(z)$ is the Legendre polynomial of order $2n$. As we assume the distribution function to be a symmetric function of $p_{z}$, only polynomials of even orders are necessary. The zeroth order moment is nothing but the energy density, $\mathcal{L}_{0}=\mathcal{E}$, while the first moment gives the difference between the longitudinal and the transverse pressures, $\mathcal{L}_{1}=\mathcal{P}_{L}-\mathcal{P}_{T}$.

If one employs the relaxation time approximation $C[f]=(f-f_{\rm eq})/\tau_{R}$ for the collision term one can derive a set of coupled equations for the moments ${\cal L}_{n}$, which can be solved by truncations [6]. However, for a more realistic collision kernel, moments other than the $\mathcal{L}_{n}$’s appear, which complicates the analysis. This is the case in particular when one uses the collision term describing gluon elastic collisions in the small-angle approximation [9, 11]. As shown in Appendix A, the equations for the moments $\mathcal{L}_{n}$ are then coupled to other kinds of moments, and no obvious and simple scheme emerges that would allow us to close the resulting system of equations.

Of course, the Boltzmann equation (1) can be solved directly, with the elastic collision term treated within the small-angle approximation (see e.g. [9, 14]), or beyond [2, 15], as well as with both elastic and inelastic processes included in the collision kernel (see e.g. Refs. [1, 3, 16]). However these approaches are very demanding numerically. In this paper, we then attempt an intermediate approach between the moment method of Ref. [6] and the direct analysis of the full Boltzmann equation. We make a mode expansion with respect to the angular variable $\cos\theta=p_{z}/p$ while the dependence on the absolute value of the 3-dimensional momentum $p=|{\boldsymbol{p}}|$ is treated fully. As we shall see, since in most cases of practical interest a few mode functions are sufficient to capture the main physics, this represents a substantial gain: a two-dimensional kinetic equation is replaced by a finite (small) set of coupled one-dimensional equations.111To appreciate the numerical gain, note that in Ref. [14], the angular variable is discretized on a grid with 256 points, whereas only a few modes functions are needed to get an accurate solution. Aside from this technical aspect, the mode functions may represent an interesting set of degrees of freedom that provide useful insight into the dynamics of the expanding system.

We consider the following expansion of the momentum distribution function,

[TABLE]

The mode function $f_{n}(\tau,p)$ can be extracted from the full distribution function by

[TABLE]

In terms of the mode functions, the moments (3) are given by

[TABLE]

The particle number density ${\cal N}$ and the energy density $\mathcal{E}$ depend only on the zeroth mode function $f_{0}(\tau,p)$, and are given respectively by

[TABLE]

and

[TABLE]

where $\nu_{\text{g}}=2(N_{c}^{2}-1)=16$ denotes the number of internal degrees of freedom of the gluons. The longitudinal and the transverse pressures depend on both the zeroth and first mode functions. In terms of the moments ${\cal L}_{0}$ and ${\cal L}_{1}$ they read (see Eq. (6))

[TABLE]

In the small-angle approximation, the collision term for the gluon-gluon elastic scattering takes the form [9, 11]

[TABLE]

where $\alpha_{s}$ is the strong coupling constant and $l_{\rm Cb}$ is the Coulomb logarithm, treated here as a constant. The integrals $I_{a}$ and $I_{b}$ are defined by

[TABLE]

and

[TABLE]

Since the integrand of $I_{b}$ is isotropic in ${\boldsymbol{p}}$ and linear in $f$, $I_{b}$ can be expressed in terms of $f_{0}$ alone:

[TABLE]

In contrast, all the mode functions contribute to the integral $I_{a}$ as it is nonlinear in $f$,

[TABLE]

By applying the mode projection operator $\frac{4n+1}{2}\int_{-1}^{1}d(\cos\theta)\,P_{2n}(\cos\theta)$ to the Boltzmann equation (1), we can derive the kinetic equations for the mode function $f_{n}(\tau,p)$

[TABLE]

where the collision term reads

[TABLE]

and we have rescaled the time variable, setting

[TABLE]

A further rescaling will be performed later in order to measure time in units of a suitably defined relaxation time (see Eq. (29) below). Here and in the following, the time argument of $f$ is often omitted to alleviate the notation. The numerical coefficients that appear in the drift terms are given by

[TABLE]

and

[TABLE]

The coefficient $A_{n,m,l}$ in the collision term is defined by

[TABLE]

The derivation of these equations and the explicit form of $A_{n,m,l}$ can be found in Appendix B.

The kinetic equations for the mode functions, Eq. (15), are coupled to each other through the drift terms ($\propto 1/\tau$), the nonlinear term in the collision kernel, as well as the integral $I_{a}$. They constitute an infinite hierarchy of equations, and in order to solve them in practice we need to implement some truncation. In Ref. [6], it is shown that such truncation works efficiently already at a low order in the moment method. In Sec. 4, we will demonstrate that a low order truncation works successfully also for the angular mode expansion. The truncation that we shall consider amounts to solving the equations for the first $N+1$ mode functions (15), i.e., for $n=0,\cdots,N$ and set to zero all the mode functions with $n>N$.

The conservation laws are not affected by the truncation, except the energy conservation law in the special case of the truncation at order $N=0$. These conservation laws are easily deduced by integrating the kinetic equation with suitable weights. The conservation law of the number density reads

[TABLE]

It involves only the $f_{0}$ mode function, and its form is not modified in any truncation. In contrast, the conservation law of the local energy density

[TABLE]

involves also the $f_{1}$ mode through the longitudinal pressure $\mathcal{P}_{L}$. In the truncation orders $N\geq 1$, its form is not affected. However, in the truncation order $N=0$, the first order moment (i.e. the pressure anisotropy) is forced to vanish, implying that the energy density and the pressure are always related as $\mathcal{P}_{L}=\mathcal{P}_{T}=\mathcal{E}/3$ (for massless particles). In this case, the conservation law takes the form

[TABLE]

This is nothing but the ideal hydrodynamic equation, whose solution is $\mathcal{E}(\tau)=\mathcal{E}(\tau_{0})\left(\tau_{0}/\tau\right)^{4/3}$. Of course, this is also the regime expected when all higher mode functions $n=1,2,\cdots$ have been damped by the collisions, and where the evolution of the system is fully described by the zeroth mode function.

3 Initial conditions and time scales

Before we present some numerical results, we discuss shortly the initial conditions used in the calculations. We also provide an estimate of the time scale that controls the collision rate.

3.1 Initial conditions

In a collision of high-energy nuclei, gluons whose typical momentum scale is given by the saturation momentum $Q_{s}$ are emitted in the central rapidity region [17]. As in previous related studies (see e.g. [8]), we model the distribution of these gluons by the following initial distribution,

[TABLE]

where $\Theta(z)$ is the Heaviside step function, $f_{\rm ini}$ is a parameter that characterizes the initial occupation number, and $\xi_{0}>1$ measures the initial anisotropy. This particular form allows for a simple determination of the initial mode functions. We emphasize that none of our main conclusions depends crucially on the shape of this initial distribution. The time scale for the formation of the initial gluons is $\tau_{0}\sim Q_{s}^{-1}$, and throughout we shall take this time as the starting time for the kinetic evolution. Note that $Q_{s}$ is the only momentum scale, and we shall consistently express all momenta in units of $Q_{s}$.

The initial conditions for the mode functions $f_{n}(p)$ that correspond to the initial distribution (24) are given by

[TABLE]

These are plotted in Fig. 1 for several values of $n$ and an anisotropy parameter $\xi_{0}=1.5$. The mode functions for $n\geq 1$ are nonzero only in the region $Q_{s}/\xi_{0}<p<Q_{s}$, where the distribution shows anisotropy. They exhibit oscillations whose wave number grows with increasing $n$.

In Fig. 2, we show how the original distribution is reconstructed by summing up the mode functions for several truncation orders. More precisely, we plot $\sum_{n=0}^{N}f_{n}(t,p)P_{2n}(\cos\theta)$ at $p_{z}=0$ ($\cos\theta=0$) for several $N$. For large $N$, this sum indeed approaches the original step function. Not surprisingly, the convergence is slower at the corner of the distribution, near the sharp discontinuity at $p=Q_{s}$. However, we shall see that such small structures in the distribution functions are quickly smeared out by collisions.

3.2 Time scales

To analyze the results of the numerical calculations to be presented in the next section, it is useful to measure the time in units of a quantity that can be associated to a basic relaxation time $\tau_{R}$. This time can be roughly estimated by comparing the collision term (10) with a relaxation time approximation, $(f-f_{\rm eq})/\tau_{R}$. From this comparison, we find

[TABLE]

The combination $4\pi N_{c}^{2}\alpha_{s}^{2}\,l_{\rm Cb}$ is that already encountered in the definition of $\tilde{\tau}$ in Eq. (17). Thus defined, the relaxation time is a momentum and also time-dependent quantity. We fix the momentum dependence by evaluating $\tau_{R}$ at the typical momentum scale of the gluons at the early time, i.e. $p=Q_{s}$. The time dependence comes from the integral $I_{a}$ defined in Eq. (11). If the typical occupation number is not much larger than unity, one can neglect the Bose-enhancement factor.222We shall not consider in this paper the case of large occupation, $f_{\rm ini}\sim 1/\alpha_{s}$, which has been thoroughly analyzed in Ref. [14]. (We typically use for initial occupation number $f_{\rm ini}=0.1$.) In that case, the integral $I_{a}$ is approximately equal to the number density ${\cal N}$ (modulo the factor $\nu_{\text{g}}$) and, for the expanding case where ${\cal N}(\tau)\propto 1/\tau$, it can be evaluated as

[TABLE]

The relaxation time for the expanding case is therefore estimated as

[TABLE]

Since the system gets diluted proportionally to $\tau$ by the longitudinal expansion, the expansion rate decreases as $1/\tau$. With the present estimate of $\tau_{R}(\tau)$, the ratio between the collision rate and the expansion rate is therefore independent of $\tau$. This is an oversimplification however. In the following we shall present results in terms of a time variable measured in units of the relaxation time $\tau_{R}$ defined by the expression (28) at the initial time $\tau_{0}$, that is

[TABLE]

By doing so, we do not need to specify the values of the various parameters entering the definition (29), in particular the coupling constant $\alpha_{s}$. Note however that for typical values $\alpha_{s}\simeq 0.3$, $\xi_{0}\simeq 1$, $l_{\rm Cb}\simeq 1$, one finds $\tau_{0}/\tau_{R}\simeq 0.2f_{\rm ini}$, as a crude order of magnitude.

In the next section, we shall also consider the non expanding case. In this case, the density is constant, and so is $I_{a}$. The relaxation time is then independent of time and is defined as in Eq. (29), with $\tau_{0}$ replaced by $1/Q_{s}$,

[TABLE]

4 Numerical results

In this section, we show numerical results obtained by solving the kinetic equations (15) for the mode functions with the initial conditions (25). To show the efficiency of the angular mode expansion method, we compare these results with those obtained by solving the full Boltzmann equation, Eqs. (1) and (10), using the method presented in Ref. [14]. Since the kinetic equations for the individual mode functions take almost the same form as a Boltzmann equation for an isotropic system, they can be solved by the same method as the original Boltzmann equation but with much lower numerical cost, since they are one-dimensional instead of two-dimensional.

Below we employ a fixed set of parameters for the initial conditions $f_{\rm ini}=0.1$ and $\xi_{0}=1.5$. In numerical calculations, we have nonuniformly discretized the momentum variable $p\in[p_{\text{min}},p_{\text{max}}]$ into 200–400 grid points, where $p_{\text{max}}/Q_{s}=4$–8 and $p_{\text{min}}/Q_{s}=10^{-2}$. The time evolution is solved by an implicit method.

We shall first consider the non expanding case and analyze how the system isotropize and thermalize. Then we turn to the expanding case, and consider the same issues for various values of the collision rate.

4.1 Non-expanding plasma

For a non-expanding plasma, we need to solve the following set of kinetic equations

[TABLE]

where the drift terms ($\propto 1/\tau$) in Eq. (15) have been left aside. We shall do so by truncating this infinite set of equations to a finite number, $n\leq N$. We shall see that values $N\leq 4$ are actually sufficient in order to get accurate results (in most cases, even $N=1$ is enough). The results obtained with this method will be compared with the solution of the corresponding kinetic equation solved by the method of Ref. [14]. In the non-expanding system, the value of the initial time does not affect the results, so we chose $t_{0}=0$. We shall express all results as a function of $t/\tau_{R}$, where $\tau_{R}$ is the relaxation time defined in Eq. (30).

Both particle number and energy are conserved to an accuracy better than $10^{-3}$ within the time range shown in the following figures.

For a first orientation, we show in Fig. 3, left panel, the longitudinal and the transverse pressures divided by the energy density as a function of time. One sees that the pressures are isotropized on a time scale $\lesssim\tau_{R}$. The right panel of Fig. 3 shows the time evolution of the moments $\mathcal{L}_{n}$ defined in Eq. (3). As was observed already in Ref. [8], there is a hierarchy of the moments $|\mathcal{L}_{0}|>|\mathcal{L}_{1}|>|\mathcal{L}_{2}|>\cdots$, and all these moments decay quickly in time, the higher the rank $n$ the faster the decay. This is another and more detailed way to visualize the isotropization of the system as a function of time. The results in Fig. 3 were obtained by solving the full Boltzmann equation. However as we shall see now, equivalent results are obtained with low order truncations.

In the non-expanding system, the different angular modes couple to each other only through the terms nonlinear in $f$ that are present in the collision kernel of Eq. (16). One of such nonlinearity is hidden in the integral $I_{a}$ given in Eq. (14). The time evolution of $I_{a}$ (divided by the number density $\cal N$) is plotted for several truncation orders $N\leq 4$ in Fig. 4, left panel. For comparison, the result obtained from the full Boltzmann equation is also plotted as a dotted line. The nonlinearity is measured here by the deviation of $I_{a}/{\cal N}$ from unity (when $f\ll 1$, $I_{a}\simeq{\cal N}$). Small deviations from the full result are seen only at early times and rapidly disappear in higher order truncations. At later time, even the 0th order truncation agrees with the full result. The right panel of Fig. 4 displays the time evolution of the pressure ratio $\mathcal{P}_{L}/\mathcal{P}_{T}$. Note that in the lowest truncation order $N=0$, the anisotropy is entirely neglected and the pressure ratio is always one. Truncations with $N\geq 1$ fully account for the anisotropy of pressures. As seen in the figure, all the results for $N=1,2,4$ nicely agree with the full result. This indicates that the coupling to higher modes $n=2,3,\cdots$ does not play much of a role in the time evolution of the pressures: it is enough to keep the lowest angular modes $f_{0}$ and $f_{1}$ to get an accurate description.

To have a closer look at the efficiency of the truncation, we plot in Fig. 5 the first two mode functions $f_{0}(p)$ and $f_{1}(p)$ as a function of momentum at some early time $t/\tau_{R}=0.01$, for several truncation orders. The results corresponding to different truncations completely overlap with each other, and are in perfect agreement with the full result represented by the dotted line. Thus, already at early time ($t\ll\tau_{R}$) the truncation of the angular modes works extremely well. Note that this does not necessarily imply that all the other angular modes are vanishingly small. Simply they do not significantly affect the dynamics of the lowest modes, in which we are primarily interested.

We can push this analysis a bit further and make a more detailed comparison of the evolution of the mode functions $f_{n}(p)$ for different $n$. In Fig. 6, the mode functions for $n=0,1,2,3,4$ at times $t/\tau_{R}=0$, 0.01, 0.1, 0.3 are plotted as a function of momentum. These results are obtained with the truncation $N=4$. As we have already seen in Sec. 3.1, the angular modes $n\geq 1$ exhibit initially rapid oscillations. However, these oscillations are quickly washed out, as seen in the right upper panel ($t/\tau_{R}=0.01$). As time goes on, all the higher mode functions decay in the whole momentum region and only the mode $f_{0}$ survives at a later time. In the lower two panels, the thermal distribution $f_{\rm eq}(p)=1/(e^{(p-\mu)/T}-1)$ is also plotted as a gray dashed line. In the non-expanding system, with conserved particle number, one can calculate the temperature $T$ and chemical potential $\mu$ of the final state from the initial condition. For our initial condition, these are $T=0.23\,Q_{s}$ and $\mu=-0.075\,Q_{s}$. From $t/\tau_{R}\simeq 0.3$, when all the higher mode functions have essentially decayed, the zero mode function nicely agrees with the thermal distribution.

The thermalization of the system involves not only the decay of the modes $f_{n}(p)$ for $n\geq 1$, a process that we may refer to as isotropization, but also the rearrangement of the various isotropic momentum modes, that is the evolution of $f_{0}(p)$ towards a thermal distribution, as we have just mentioned333Note that we do not find particular distinct time scales that could be attributed to each of these processes which rather progress in parallel.. To study more precisely this phenomenon, we have extracted an effective temperature $T_{\text{eff}}$ and chemical potential $\mu_{\text{eff}}$ by fitting the zero mode function $f_{0}(p)$ to an equilibrium distribution $f_{\rm eq}(p)=1/(e^{(p-\mu_{\text{eff}})/T_{\text{eff}}}-1)$ in the momentum range $0<p<Q_{s}$ (see e.g. [11] where a similar strategy has been used). The resulting values are plotted in Fig. 7 as a function of time. Except at early times t/\tau_{R}\raisebox{-2.15277pt}{>\stackrel{{\scriptstyle\textstyle<}}{{\sim}}>}0.1, the fit is accurate, thus we do not indicate fitting errors in the figure. Both the effective temperature and the chemical potential converge toward their expected values that can be deduced from the conservation laws, and which are indicated by black dashed lines. For the effective temperature, we have employed also another definition $T_{\text{eff}}=I_{a}/I_{b}$, which coincides with the actual temperature if the distribution function is an equilibrium distribution. Both definitions of the effective temperature yield nearly the same values at late time. As can be seen in Fig. 8 the agreement of $f_{0}(p)$ with a thermal distribution is excellent already at $t/\tau_{R}\simeq 0.3$. By that time, the pressures are nearly equilibrated, but not quite, the transverse pressure exceeding the longitudinal pressure by a few percents (see Figs. 3 and 4). Part of the reason is that, by this time the high momentum modes are not fully thermalized. This is obviously the case at $t/\tau_{R}=0.1$, as Fig. 8 shows (the low momentum modes quickly thermalize).

4.2 Expanding plasma

We turn now to the expanding case. Our goal here is twofold. On the one hand, we wish to verify that the convergence of the angular mode expansion remains of a good quality. And on the other hand, we want to verify the conclusions of another study, based on the analysis of a set of coupled equations for special moments, and a simplified collision kernel based on a relaxation time approximation [6, 7]. As shown in the study in question, the competition between free streaming and collisional effects can be understood from the competition between two fixed points that govern the evolution of the moments. The first fixed point is that of the free streaming motion, and it is independent of the collision kernel. The main property of this fixed point is that it is mostly determined by the two equations that control the evolution of the lowest moments ${\cal L}_{0}$ and ${\cal L}_{1}$, that is, in the present case, it involves only the two mode functions $f_{0}$ and $f_{1}$. The other moments, or mode functions, correct the position of the fixed point, but only in a moderate way, and most importantly they do not alter the qualitative behavior of the solution. The other fixed point corresponds to the hydrodynamic behavior to be expected at late time, unless the collision rate is too small as compared to the expansion rate. The hydrodynamics is dominated by the first two mode functions, or equivalently the two moments ${\cal L}_{0}$ and ${\cal L}_{1}$ that contribute directly to the energy momentum tensor.

We shall start by discussing the case where the collision rate is comparable to the expansion rate, that is $\tau_{0}/\tau_{R}=1$. As we shall see, in this case, the collisions are nearly as efficient as in the non expanding case, they isotropize the system and drives it to local equilibrium. Then we discuss a few results obtained for smaller collision rates, where the evolution of the system becomes eventually dominated by the free streaming fixed point.

4.2.1 Large collision rate ($\tau_{0}/\tau_{R}=1$)

By large collision rate, we mean here that the collision rate is comparable to the expansion rate. The expansion plays an important role, but to a large extent the features that we have exhibited in the previous subsection for the non-expanding case survive.

The plots in Fig. 9 reveal the quality of the truncations in the expanding case. It is to be compared to the corresponding plot for the non-expanding case in Fig. 4: the quality of the truncation is in both cases the same, and the truncation $N=1$ becomes quickly accurate (note however the change of horizontal scale between the two panels, from logarithmic to linear).

The accuracy of the truncation is further confirmed by the analysis of the mode functions. As an example, we display in Fig. 10 the mode functions $f_{0}$ and $f_{1}$ as a function of momentum, at time $\tau-\tau_{0}=0.5\,\tau_{R}$. Clearly, even the truncation $N=1$ reproduces accurately the values of these functions over the entire momentum range. As we did in the non expanding case, we can define an effective temperature $T_{\rm eff}$ and an effective chemical potential $\mu_{\rm eff}$ by fitting $f_{0}(p)$ to a local equilibrium distribution function. The values of the effective parameters thus obtained are plotted in Fig. 11. The time-dependence of the effective temperature is compared to that expected for a system in local equilibrium, namely $T(\tau)\sim\tau^{-1/3}$ (ideal hydrodynamical behavior). One sees that this behavior is well reproduced for times $\tau-\tau_{0}\gtrsim\tau_{R}$.

The mode function $f_{0}(p)$ is plotted (multiplied by $p^{2}$) in Fig. 12. This evolution follows a smooth pattern, and in fact is essentially that of an equilibrium distribution function with slowly evolving parameters $T_{\rm eff}$ and $\mu_{\rm eff}$. As can be seen on this figure, by time $\tau-\tau_{0}\gtrsim\tau_{R}$, $f_{0}(p)$ is accurately described by a local equilibrium distribution. However, it is important to note that the mode function $f_{1}(p)$ does not vanish at this time, and indeed the longitudinal pressures are not yet equilibrated by this time as can be seen in Fig. 9, right panel.

4.2.2 Smaller collision rates $(\tau_{0}/\tau_{R}=0.1,\ \tau_{0}/\tau_{R}=0.01)$

We now examine how the system evolves when the collision rate is small compared to the expansion rate. We shall consider two cases.

In the first case, $\tau_{0}/\tau_{R}=0.1$, the results obtained differ only slightly from those of the case $\tau_{0}/\tau_{R}=1$ that we have just examined: the truncations loose accuracy, but even the truncation $N=1$ still yields results that compare well with the exact solution. In particular, isotropization and thermalization are clearly achieved at late times. This is illustrated in the left panels of Figs. 13 and 16, which compare well with the corresponding figures in the previous subsection.

This is also confirmed by the plot of the mode functions $f_{0}$ and $f_{1}$ displayed in Fig. 14, as well as those of the effective temperature and chemical potential.

When $\tau_{0}/\tau_{R}=0.01$, the second case that we consider, we enter a regime where more mode functions are needed to reproduce the exact results. This case is illustrated by the plots in the right-hand panels of Figs. 13 and 16. The truncation $N=1$ is no longer sufficient, as indeed some artefacts appear in this case: for instance the longitudinal pressure becomes negative (see the right bottom panel of Fig. 13), an unphysical feature already observed in the two-moment truncation of Ref. [7]. Note, however, that from a practical point of view, the expansion in mode functions remains very useful since in all cases only a few modes are needed (in the present case $N=4$ is sufficient to reproduce with good accuracy all the quantities considered, see e.g. Fig. 15. We note also that in this case of small collision rate, the thermalization of the system is not achieved during the time span of the calculation. In particular the mode function $f_{0}(p)$ is never well approximated by a local equilibrium function as can be seen in the right panel of Fig. 16. It is also clear from the plot of $T_{\rm eff}$ in Fig. 15 that the hydrodynamic regime is not reached.

5 Conclusions

We have shown that the expansion of the distribution function in angular modes allows for a significant reduction of the numerical effort involved in solving the Boltzmann equation for longitudinally expanding quark-gluon plasmas. This is because only a few angular modes are enough to capture the main features of the time evolution of the energy density and the pressures, which are related to the lowest two moments of the distribution function, i.e. to the first two angular modes $f_{0}$ and $f_{1}$. The number of modes needed for an accurate description of the distribution function depends on the ratio between the expansion rate and the collision rate. For most cases of physical interest, the first two modes suffice, higher modes being needed only in cases of small collision rates, i.e., when the evolution of the system is strongly dominated by free streaming.

We have presented results of numerical calculations, first for a non expanding system, and then for a boost invariant, longitudinally expanding system. In both cases we found a rapid convergence of the expansion. In the non-expanding case, we observe the isotropization of the distribution, measured by the decay of the various mode functions $f_{n\geq 1}$. The thermalization is observed though the fact that the spherical mode $f_{0}(p)$ is well described at late time by a thermal distribution, with slowly evolving temperature and chemical potential. In the expanding system, the expansion delays the effect of the collisions. We find however that this delay alters the qualitative pattern only when the collision rate is much smaller than the expansion rate. When both rates are comparable, we find that the system approaches equilibrium in very much the same way as in the non-expanding case. However, the pressure anisotropy takes more time to disappear, in particular it subsists even after the spherical mode is well described by a thermal distribution. This is due to the fact that $f_{0}$ and $f_{1}$ remain strongly coupled, and the expansion hinders the decay of $f_{1}$. When the collision rate is too small, the thermalization itself is hindered by the fast expansion.

These conclusions strengthen those obtained in [6, 7], based on specific moments of the distribution function and a relaxation time approximation for the collision kernel. The present study should be extended to more realistic collision kernel involving both elastic and inelastic processes. This is work for the future.

Acknowledgments

We thank Li Yan for useful exchanges at an early stage of this project. We also thank François Gelis for comments on the manuscript.

A Equations for the moments

In this Appendix, we derive the equation of motion for the moments $\mathcal{L}_{n}$ defined in Eq. (3) from the Boltzmann equation (1) with the small-angle approximation for the collision term (10).

From the Boltzmann equation (1), it follows that

[TABLE]

By applying integration by parts and the recursion relations for the Legendre polynomials

[TABLE]

and

[TABLE]

one can compute the moments of the drift term as [6]

[TABLE]

with coefficients

[TABLE]

These coefficients are related with those given by Eqs. (18-19) as $\alpha_{n}=4a_{n}-\tilde{a}_{n}$, $\beta_{n}=4b_{n}-\tilde{b}_{n}$ and $\gamma_{n}=4c_{n}-\tilde{c}_{n}$.

The moment of the collision term can be computed in a similar way. We find

[TABLE]

The momentum integrals that appear in the right hand side are different kinds of moments from $\mathcal{L}_{n}$. For those, let us introduce

[TABLE]

and

[TABLE]

Finally, the equation for the moments $\mathcal{L}_{n}$ reads

[TABLE]

This equation is clearly not closed within the moments $\mathcal{L}_{n}$. If one tries deriving similar equations for the moments $\mathcal{J}_{n}$ and $\mathcal{K}_{n}$, one encounters new kinds of moments. This makes it difficult to compute the dynamics of this system based on simple equations for a closed set of moments.

B Derivation of the kinetic equations for the angular mode functions

By applying the mode projection operator $\frac{4n+1}{2}\int_{-1}^{1}d(\cos\theta)\,P_{2n}(\cos\theta)$ to the Boltzmann equation (1), we can derive the equations of motion for the mode function $f_{n}(p)$. First, we rewrite the Boltzmann equation in terms of the variables $p=|{\boldsymbol{p}}|$ and $x=\cos\theta$,

[TABLE]

with

[TABLE]

From Eqs. (A2), (A3) and the orthogonality relation

[TABLE]

one can derive the following formulas:

[TABLE]

and

[TABLE]

By using these formulas, one can find the projection of the expansion drift term as is given in Eq. (15).

Next, we compute the projection of the collision term. The term proportional to $I_{a}$ is simple because the Legendre polynomials are eigen functions of the Laplacian.

[TABLE]

The term proportional to $I_{b}$ is much more involved since it is nonlinear in $f$.

[TABLE]

where we have defined

[TABLE]

The explicit form of $A_{n,m,l}$ can be derived by using the formula [Eq. 8.915.5 (p.986) of Ref. [18]]

[TABLE]

where

[TABLE]

Here, $\alpha_{0}$ should be understood as 1. Thanks to this formula, we can show

[TABLE]

Collecting all the above results, we arrive in Eq. (15).