{'fraction_non_alphanumeric': 0.03718106153636239, 'fraction_numerical': 0.012962204939282303, 'mean_word_length': 4.43908472479901, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving the Boltzmann equation in N log N

In [C. Mouhot and L. Pareschi,"Fast algorithms for computing the Boltzmann collision operator,"Math. Comp., to appear; C. Mouhot and L. Pareschi, C. R. Math. Acad. Sci. Paris, 339 (2004), pp. 71-76], fast deterministic algorithms based on spectral methods were derived for the Boltzmann collision operator for a class of interactions including the hard spheres model in dimension three. These algorithms are implemented for the solution of the Boltzmann equation in two and three dimension, first for homogeneous solutions, then for general non homogeneous solutions. The results are compared to explicit solutions, when available, and to Monte-Carlo methods. In particular, the computational cost and accuracy are compared to those of Monte-Carlo methods as well as to those of previous spectral methods. Finally, for inhomogeneous solutions, we take advantage of the great computational efficiency of the method to show an oscillation phenomenon of the entropy functional in the trend to equilibrium, which was suggested in the work [L. Desvillettes and C. Villani, Invent. Math., 159 (2005), pp. 245-316].


Introduction
The construction of approximate methods of solution for the Boltzmann equation has a long history tracing back to D. Hilbert, S. Chapmann and D. Enskog [13] at the beginning of the last century. The mathematical difficulties related to the Boltzmann equation make it extremely difficult, if not impossible, the determination of analytic solutions in most physically relevant situations. Only in recent years, starting in the 70s with the pioneering works by A. Chorin [14] and G. Sod [46], the problem has been tackled numerically with particular care to accuracy and computational cost. Even nowadays the deterministic numerical solution of the Boltzmann equation still represents a challenge for scientific computing.
Most of the difficulties are due to the multidimensional structure of the collisional integral, since the integration runs on a highly-dimensional unflat manifold. In addition the numerical integration requires great care since the collision integral is at the basis of the macroscopic properties of the equation. Further difficulties are represented by the presence of stiffness, like the case of small mean free path [24] or the case of large velocities [19].
For such reasons realistic numerical simulations are based on Monte-Carlo techniques. The most famous examples are the Direct Simulation Monte-Carlo (DSMC) methods by Bird [3] and by Nanbu [35]. These methods guarantee efficiency and preservation of the main physical properties. However, avoiding statistical fluctuations in the results becomes extremely expensive in presence of non-stationary flows or close to continuum regimes. Among deterministic approximations, perhaps the most popular method is represented by the so-called Discrete Velocity Models (DVM) of the Boltzmann equation. These methods [28,45,6,10,38] are based on a cartesian grid in velocity and on a discrete collision mechanism on the points of the grid that preserves the main physical properties. Unfortunately DVM are not competitive with Monte-Carlo methods in terms of computational cost and their accuracy seems to be less than first order [36,37,38]. In this work we are interested in high-order deterministic methods and therefore we shall not discuss algorithms based on DVM, and we refer the reader to the work in preparation [33].
More recently a new class of numerical methods based on the use of spectral techniques in the velocity space has been developed. The methods were first derived in [40], inspired from spectral methods in fluid mechanics [11] and by previous works on the use of Fourier transform techniques for the Boltzmann equation (see [5] for instance). The numerical method is based on approximating the distribution function by a periodic function in the phase space, and on its representation by Fourier series. The resulting Fourier-Galerkin approximation can be evaluated with a computational cost of O(n 2 ) (where n is the total number of discretization parameters in velocity), which is lower than that of previous deterministic methods (but still larger then that of Monte-Carlo methods).
It was further developed in [41,43] where evolution equations for the Fourier modes were explicitly derived and spectral accuracy of the method has been proven. Strictly speaking these methods are not conservative, since they preserve mass, whereas momentum and energy are approximated with spectral accuracy. This trade off between accuracy and conservations seems to be an unavoidable compromise in the development of numerical schemes for the Boltzmann equation (with the noticeably exception of [39]) . We recall here that the spectral method has been applied also to non homogeneous situations [20,22], to the Landau equation [19,42], where fast algorithms can be readily derived, and to the case of granular gases [34,21]. For a recent introduction to numerical methods for the Boltzmann equation and related kinetic equations we refer the reader to [15]. Finally let us mention that A. Bobylev & S. Rjasanow [7,8] have also constructed fast algorithms based on a Fourier transform approximation of the distribution function.
In [32,31] a fast spectral method was proposed for a class of particle interactions including pseudo-Maxwell molecules (in dimension 2) and hard spheres (in dimension 3), on the basis of the previous spectral method together with a suitable semi-discretization of the collision operator. This method permits to reduce the computational cost from O(n 2 ) to O(n log 2 n) without loosing the spectral accuracy, thus making the method competitive with Monte-Carlo. The principles and basic features of this new method will be presented in the next sections.
The rest of the paper is organized as follows. Section 2 is devoted to a short introduction on the Boltzmann equation and its physical properties. Next in Section 3 we explain the principles of the different spectral algorithms used to compute the collision operator. Several numerical results and comparisons to exact solutions as well as to Monte-Carlo methods are given in Section 4. An application to a challenging non homogeneous test case is finally given in Section 5. Some final considerations close the paper in the last Section.

The Boltzmann equation
The Boltzmann equation describes the behavior of a dilute gas of particles when the only interactions taken into account are binary elastic collisions. It reads for is the time-dependent particles distribution function in the phase space. The Boltzmann collision operator Q is a quadratic operator local in (t, x). The time and position acts only as parameters in Q and therefore will be omitted in its description The collision kernel B is a non-negative function which by physical arguments of invariance only depends on |v − v * | and cos θ =ĝ In this work we are concerned with short-range interaction models. More precisely we assume that B is locally integrable. This assumption is satisfied by the hard spheres model, which writes in dimension d = 3 and is known as Grad's angular cutoff assumption when it is (artificially) extended to interactions deriving from a power-law potentials. As an important benchmark model for the numerical simulation we therefore introduce the so-called variable hard spheres model (VHS), which writes for some γ ∈ [0, 1] and a constant C γ > 0.
For this class of model, one can split the collision operator as Boltzmann's collision operator has the fundamental properties of conserving mass, momentum and energy The functional − f log f is the entropy of the solution. Boltzmann's H theorem implies that any equilibrium distribution function, i.e., any function which is a maximum of the entropy, has the form of a locally Maxwellian distribution where ρ, u, T are the density, mean velocity and temperature of the gas, defined by For further details on the physical background and derivation of the Boltzmann equation we refer to [13] and [49].

The spectral methods
In this section we shall explain the principles of the algorithms to compute the collision integral for a fixed value of the spatial variable x. Indeed it is well-known that one can reduce to this case by some splitting strategy (see [41,20] for example).

3.1.
A general framework. We consider the spatially homogeneous Boltzmann equation written in the following general form where Q is given by . In the equations above, C is some given unbounded domain, and Θ, Θ ′ , Θ ′ * are suitable functions, to be defined later. This general framework emphasizes the translation invariance property of the collision operator, which is crucial for the spectral methods. We will be more precise in the next paragraphs for some changes of variables allowing to reduce the classical operator (2.1) to the form (3.2).
A problem associated with deterministic methods which use a fixed discretization in the velocity domain is that the velocity space is approximated by a finite region. Physically the domain for the velocity is R d , and the property of having compact support is not preserved by the collision operator. In general the collision process spreads the support by a factor √ 2 in the elastic case (see [44,30] and also [29] for similar properties in the inelastic case). As a consequence, for the continuous equation in time, the function f is immediately positive in the whole domain R d . Thus, at the numerical level, some non physical condition has to be imposed to keep the support of the function in velocity uniformly bounded. In order to do this there are two main strategies.
• One can remove the physical binary collisions that will lead outside the bounded velocity domain, which means a possible increase of the number of local invariants (that is the functions ϕ such that (ϕ ′ * + ϕ ′ − ϕ * − ϕ) is zero everywhere on the domain). If this is done properly (i.e., without removing too many collisions), the scheme remains conservative (and without spurious invariants). However, this truncation breaks down the convolution-like structure of the collision operator, which requires the translation invariance in velocity. Indeed the modified collision kernel depends on v through the boundary conditions. This truncation is the starting point of most schemes based on Discrete Velocity Models in a bounded domain.
• One can add some non physical binary collisions by periodizing the function and the collision operator. This implies the loss of some local invariants (some non physical collisions are added). Thus the scheme is not conservative anymore, except for the mass if the periodization is done carefully. In this way the structural properties of the collision operator are maintained and thus they can be exploited to derive fast algorithms. This periodization is the basis of the spectral methods.
Therefore, we consider the space homogeneous Boltzmann equation in a bounded domain in velocity D T = [T ; T ] d (0 < T < ∞). We need to truncate the integration in y and z since periodization would yield infinite result if not. Thus we set y and z to belong to some truncated domain C R ⊂ C (the parameter R refers to its size and will be defined later). For a compactly supported function with support included in B S , the ball centered at 0 with radius S > 0, one has to prescribe suitable relations (depending on the precise change of variable and truncation chosen) between S, R, T in order to retain all possible collisions and at the same time prevent intersections of the regions where f is different from zero (dealiasing condition). Then the truncated collision operator reads for v ∈ D T (the expression for v ∈ R d is deduced by periodization). By making some changes of variable on v, one can easily prove for the two choices of variables y, z of the next subsections, that for any function ϕ periodic on D T the following weak form is satisfied Now, we use the representation Q R to derive spectral methods. Hereafter, we use just one index to denote the d-dimensional sums with respect to the vector k = (k 1 , .., k d ) ∈ Z d , hence we set The approximate function f N is represented as the truncated Fourier series In a Fourier-Galerkin method the fundamental unknowns are the coefficientsf k , k = −N, . . . , N. We obtain a set of ODEs for the coefficientsf k by requiring that the residual of (3.3) be orthogonal to all trigonometric polynomials of degree less than N. Hence for k = −N, . . . , N By substituting expression (3.5) in (3.4) we get where (3.9) β(l, m) = C R B(y, z)e i π T l·Θ ′ (y,z)+m·Θ ′ * (y,z) dy dz.
The spectral equation is the projection of the collision equation in IP N , the (2N + 1) d -dimensional vector space of trigonometric polynomials of degree at most N in each direction, i.e., where P N denotes the orthogonal projection on IP N in L 2 (D T ). A straightforward computation leads to the following set of ordinary differential equations on the Fourier coefficients with the initial condition

3.2.
Classical spectral methods. In the classical spectral method [41], a simple change of variables in (2.1) permits to write Finally B c is defined by The Boltzmann operator (3.12) is now written in the form (3.2) with (y, z) = (g, ω) ∈ R d ×S d−1 =: C. Moreover, from the conservation of the momentum v ′ * +v ′ = v * + v and the energy |v ′ * | 2 + |v ′ | 2 = |v * | 2 + |v| 2 , we get the following result [40], assuming supp f ⊂ B S , • the collision operator is then given by As a consequence of this result, in order to write a spectral approximation to and extend it by periodicity to a periodic function on [−T, T ] d . We truncate the domain for (y, z) = (g, ω) as C R = B R × S d−1 for R > 0 (defining thus Q R ). Following the previous discussion on the dealiasing condition, we take R = 2S and the shortest period can be restricted to [−T, T ] d , with T ≥ (3 + √ 2)S/2 (see for a more detailed discussion [41]).
Then, we apply the spectral algorithm (3.7) and (3.8) and get the following kernel modes β c (l, m) We refer to [41,22] for the explicit computation of Fourier coefficients β c (l, m) in the VHS case where B is given by (2.3). Now, the evaluation of the right-hand side of (3.10) requires exactly O(N 2d ) operations. We emphasize that the usual cost for a DVM method based on where M is the numbers of angle discretizations.

Fast Spectral methods (FSM).
Here we shall approximate the collision operator starting from a representation which conserves more symmetries of the collision operator when one truncates it in a bounded domain. This representation was used in [7,26] to derive finite differences schemes and it is close to the classical Carleman representation (cf. [12]). The basic identity we shall need is (for u ∈ R d ) Using (3.16) the collision operator (2.1) can be written as Thus, the collision operator is now written in the form (3.2) with (y, z) ∈ R d ×R d =: for the distribution f , and the bounded domain B R × B R for (y, z) (for some R > 0). If f has support included in B S , S > 0, geometrical arguments similar to the one for the classical spectral methods (see [41,32,31]) show that we can take R = √ 2S and T as in the classical spectral method to get all collisions and prevent intersections of the regions where f is different from zero. The (truncated) operator now reads (3.18) for v ∈ D T . This representation of the collision kernel yields better decoupling properties between the arguments of the operator. From now, we can apply the spectral algorithm (3.7) and (3.8) to this collision operator and the corresponding kernel modes are given by In the sequel we shall focus on β f , and one easily checks that β f (l, m) depends only on |l|, |m| and |l · m|.
Remark 3.1. Note that the classical spectral method originates the following form of the kernel modes in the y, z notation One can notice that the condition |y + z| 2 = |y| 2 + |z| 2 ≤ R 2 couples the modulus of y and z, such that the ball is not completely covered (for instance, if y and z are orthogonal both with modulus R, the condition is not satisfied, since |y + z| = √ 2R). This explains the better decoupling properties between the argument of the collision operator of this representation.

Fast algorithms
The search for fast deterministic algorithms for the collision operator, i.e., algorithms with a cost lower than O(N 2d+ǫ ) (with typically ǫ = 1 for DVM, or ǫ = 0 for the classical spectral method), consists mainly in identifying some convolution structure in the operator (see for example [8,42]). The aim is to approximate each β f (l, m) by a sum This gives a sum of A discrete convolutions and so the algorithm can be computed in O(A N d log 2 N) operations by means of standard FFT techniques [11]. To this purpose we shall use a further approximated collision operator where the number of possible directions of collision is reduced to a finite set. We start from representation (3.18) and write y and z in spherical coordinates (note that thanks to the orthogonality condition imposed by the Dirac mass on y and z, B f depends only on the modulus of y and z). Let us denote by A a discrete set of orthogonal couples of unit vectors (e, e ′ ), which is even, i.e., (e, e ′ ) ∈ A implies that (−e, e ′ ), (e, −e ′ ) and (−e, −e ′ ) belong to A (this property on the set A is required to preserve the conservation properties of the operator). Now we define Q R,A to be where dA denotes a discrete measure on A which is also even in the sense that dA(e, e ′ ) = dA(−e, e ′ ) = dA(e, −e ′ ) = dA(−e, −e ′ ). It is easy to check that Q R,A has the same conservation properties as Q R . We make the decoupling assumption that This assumption is obviously satisfied if B f is constant. This is the case of Maxwellian molecules in dimension d = 2, and hard spheres in dimension d = 3 (the most relevant kernel for applications). Extensions to more general interactions are discussed in [32]. Let us describe the method in dimension d = 3 with B f satisfying the decoupling assumption (4.20) (see [32] for other dimensions). First we change to spherical coordinates and then we integrate first e ′ on the intersection of the unit sphere with the plane Thus we get the following decoupling formula with two degrees of freedom where S 2 + denotes the half-sphere and and Π e ⊥ is the orthogonal projection on the plane e ⊥ . In the particular case where B f = 1 (hard spheres model), we can compute the functions φ 3 R and ψ 3 Taking a spherical parametrization (θ, ϕ) of e ∈ S 2 + and taking for the set A uniform grids of respective size M 1 and M 2 for θ and ϕ we get We consider this expansion with M = M 1 = M 2 to avoid anisotropy in the computational grid. By using the Fast Fourier Transform (FFT) algorithm, the computational cost of the algorithm is then of a direct discretization on the grid for a DVM method, and O(N 6 ) of the classical spectral method).
Let us finally mention that the mathematical analysis of the fast algorithm in [32] provides the following results: • it is spectrally accurate according to the parameters N and M; • the error on the conservation laws of momentum and energy is spectrally small according to the parameter N, and no additional error (according to the speed-up parameter M) is made. This two properties were the main motivation for the development of the method of [32] described above to obtain the decomposition (4.19). Other advantages are that this particular decomposition does not introduce instability in the equation (see [32,Theorem 3.1] for instance) and it is naturally adaptative (as it is based on the rectangular quadrature rule for approximating integrals of periodic functions). Finally, another advantage of the proposed method is that it is still easy to implement since it is only based on FFT.

Numerical results in the homogeneous case
In this section we will present several numerical results for the space homogeneous equation which show the improvement of the fast spectral algorithms with respect to the classical spectral methods and how they compare with Monte-Carlo methods. The time discretization is performed by suitable Runge-Kutta methods.

Spatially homogeneous Maxwell molecules in dimension 2.
Comparison to exact solutions. We consider 2D pseudo-Maxwell molecules (i.e., the VHS model with γ = 0). In this case we have an exact solution given by with S = 1 − exp(−t/8)/2, which corresponds to the well known "BKW" solution [4]. This test is performed to check spectral accuracy, by comparing the error at a given time, when using n v = 8, 16 and 32 Fourier modes for each coordinate. We present the results obtained by the classical spectral method and the fast spectral method with different numbers of discrete angles. Figure 1 shows the relative L ∞ , L 1 , and L 2 norms of the difference between the numerical and the exact solution, as a function of time. These errors are computed according to the following formula and N = n v /2 for p = 1 and p = 2. A similar expression is used for the L ∞ error. Note that the error increases initially, and then decreases almost monotonically in time. After a long time the error starts increasing again. This effect is due to aliasing. Indeed, for a fixed computational domain, when the number of Fourier modes increases, the effect of aliasing becomes dominant over the error due to the spectral approximation. For this reason, the size of the domain is chosen in order to minimize the aliasing error. A trade-off should be obtained between aliasing and spectral error, which means that the size of the domain should be increased when increasing the number of Fourier modes. Roughly speaking, the period should be chosen in such a way that the two contributions of the error are of the same order of magnitude. In this test, the radius of the ball, which defines the computational domain is T = 4 for n v = 8, T = 5 for n v = 16 and T = 7 for n v = 32. We refer to [41] for a more detailed discussion about aliasing.
Concerning the comparison between the classical and fast spectral methods, we observe that for a fixed value of n v , the numerical error of the classical spectral method and of the fast algorithm is of the same order. Moreover, the influence of the number of discrete angles is very weak. Indeed, with only M = 4, the results are quite similar even for large n v and as expected the number of discrete angles does not affect the variations of energy, which are of the same order of magnitude as the numerical error (note that there is no variation for the momentum since in the special case of even solutions, it is preserved to 0 by the spectral scheme). In Table  1, we give a quantitative comparison of the numerical error E 1 at time T end = 1. We can also observe the spectral accuracy for the classical and fast methods: the order  Efficiency and accuracy. Now, we still consider 2D pseudo-Maxwell molecules (i.e., γ = 0) with the following initial datum where v 0 = (1, 2). In this case, we do not know the exact solution but we want to study the influence of the number of discrete angles on a non-isotropic solution.
Thus, this test is used to check the energy conservation and the evolution of highorder moments of the solution. The time step is chosen small enough to reduce the influence of the time discretization, i.e., ∆t = 0.025. Moreover, the computational domain is taken large enough with respect to the number of grid points in order to reduce the aliasing error due to the periodization of the solution. Simulations are performed with n v =16, 32 and 64 points.
In Figure 2 the relaxation of the entropy and the temperature components for the fast and classical spectral methods is shown. The energy is conserved by the continuous collision operator, but using the spectral method the total energy can change with time, it is then a good indicator on the accuracy of the numerical solution. Indeed, the total discrete energy is not exactly conserved over time, but if aliasing error is small, it is conserved within spectral accuracy, typically here the variations are about 10 −4 when n v = 32. Moreover, we observe that the number of discrete angles does not affect too much the transient regime. For instance, with only four angles on the half sphere, the relaxation of entropy and temperature components are very close to the numerical solution obtained by the classical spectral method. Finally, we plot in Figure 3 the time evolution of high-order moments of f N (t, v) given in discrete form by High-order moments give information on the accuracy of the approximate distribution function tail. Once again, we observe that the number of angles does not affect the results even if the solution is non-isotropic.
To conclude, we observe that in dimension d = 2, the fast algorithm is really efficient in terms of accuracy and computational cost compared to the classical spectral method. In Table 2 we report the computational times of the methods which show a speed-up of the fast solver independently of the number of points used in our tests and with a maximum speed-up reached for N = 64 where the fast methods with M = 4 is more than 17 times faster than the classical method.

5.2.
Spatially homogeneous hard spheres in dimension 3. In this section we consider the 3D Hard Sphere molecules (HS) model. The initial condition is chosen as the sum of two Gaussians with σ = 1 and v 0 = (2, 1, 0). The final time of the simulation is T end = 3 and corresponds approximatively to the time for which the steady state of the solution Tx with T=18, M=4, nv=32 Ty with T=18, M=4, nv=32 Tx with classic T=18, nv=32 Ty with classic T=18, nv=32      This test is used to check the evolution of moments and particularly the stress tensor P i,j , i, j = 1, · · · , 3 defined as where (u i ) i are the components of the mean velocity. As in the previous case, we compare the classical and fast methods in terms of computational time (see Table 3) and accuracy. In Figure 4, we propose the evolution of the temperature for the two methods using 32 grid points in each direction. The solution is also compared with the solution obtained from the Monte-Carlo method. The discrete temperatures agree well in this case and the efficiency of the fast algorithm is verified since the computational time is highly reduced using only M 1 = M 2 = 4 discrete angles without affecting the accuracy of the distribution function. We remark that in dimension d = 3 the speed-up of the methods becomes really evident for large values of N. Again for N = 64 and M = 4 the fast methods is more than 14 times faster.  Comparison with Monte-Carlo. Finally we compare the results obtained with our spectral method with those obtained by a Monte-Carlo scheme. We use a standard version of Monte-Carlo, which may be referred to as the Nanbu-Babowsky scheme [35,2].
In the case of Monte-Carlo methods, the moments are computed by using unbiased estimators averaged over several runs. Number of runs N runs , number of particles N p , and time step ∆t, have been chosen in such a way to balance time discretization error with statistical fluctuations.
The first run consists to take a large number of particles and to make time averaging in order to minimize the fluctuation errors: N runs = 10 3 , N p = 10 4 and ∆t = 0.01. The total computational time to compute the evolution of moments is in this case 113 min. 23 sec. On the other hand we perform a second run where we use more averaging to minimize fluctuation errors: N runs = 5. 10 3 , N p = 5. 10 3 and ∆t = 0.01. The computational time is now 100 min.
We remark that the computations have been obtained by using the hard spheres model (VHS with γ = 1), which is the most realistic. In this case, the computational time of the Monte-Carlo methods becomes larger than the case when we consider pseudo-Maxwell molecules (VHS with γ = 0) for which the collision kernel is constant and no rejection is needed. For the spectral method the computational cost is independent of the collision kernel.   From the comparison, it is obvious that, for three dimensional computations, the greater cost of the classical spectral scheme (with the same number of degrees of freedom) is compensated by a much greater accuracy, allowing better results with the same computational cost. Moreover, with the fast algorithm, the spectral scheme really becomes competitive in terms of computational time since the accuracy is not affected when we use a few number of discrete angles (for instance M = 4). In Figure 5, we compare the accuracy on the evolution of high-order moments with the different methods: the fourth order moments are very close, but the results obtained with the Monte-Carlo methods are affected by fluctuations on the tail of the distribution function, which are difficult to remove (see the evolution of the 8-th order moment). Note that for 2D and 3D pseudo-Maxwell molecules, comparisons had also been performed between Monte-Carlo methods and the classical spectral method in [41, Section 6.3].

5.3.
Stability of spectral methods with respect to non smooth data. In this subsection we perform some numerical simulations in order to study the behavior of the spectral methods when applied to non smooth data. We consider the following distribution f 0 (v) = 1 if |v| 2 ≤ 1 0 else, and use a mollified initial datum, which suitably approximates moments. We perform two numerical simulations using the fast spectral method with 32 2 and 64 2 grid points and plot the evolution of the entropy and the fourth order moment in Figure  6. Even for this discontinuous initial datum, we observe that for the two configurations the numerical entropy is decreasing and both numerical solutions converge to the same steady state. Moreover, we plot the evolution of the distribution function with respect to time in Figure 7: the fast spectral method is very stable for this numerical test even if spurious oscillations are first generated, the distribution becomes smooth and converges to an approximated Maxwellian. As expected from a Fourier-Galerkin method, the accuracy degenerates in the discontinuity region. However, surprisingly the method seems to remain stable.
Finally, let us mention that we have also performed some numerical tests when the initial datum approximates a Dirac distribution. The spectral method is still stable even if spurious oscillations are generated, but this problem is inherent to gridded methods and we refer to [23] for a rescaling method, which follows the variation of the distribution function and allows to treat concentrated distributions.

Application to the non homogeneous case
In [20,22], we performed several numerical simulations to compare the spectral scheme with Monte-Carlo methods and showed that when we are interested in the transient regime, the deterministic method becomes very efficient. Obviously, the fast spectral method will still improve the computational cost.
In the sequel, since we will be interested in the study of the trend to a global equilibrium state of the kinetic equation, we avoid the use of a splitting method by solving the whole non homogeneous equation in time by a second order Runge-Kutta method. Clearly the spectral methods apply straightforwardly to the collision operator also in this situation. The transport part is treated by the positive flux conservative method (see [18,21,20] for further details).
with periodic boundary conditions in x. We first introduce the hydrodynamical fields associated to a kinetic distribution f (t, x, v). These are the (d + 2) scalar fields of density ρ (scalar), mean velocity u (vector valued) and temperature T (scalar) defined by the formulas (2.7). Whenever f (t, x, v) is a smooth solution to the Boltzmann equation with periodic boundary conditions, one has the global conservation laws for mass, momentum and energy Therefore, without loss of generality we shall impose We shall use the following terminology: a velocity distribution of the form (6.1) will be called a Maxwellian distribution, whereas a distribution of the form will be called a local Maxwellian distribution (in the sense that the constants ρ, u and T appearing there depend on the position x). We also define the notion of relative local entropy H l , the entropy relative to the local Maxwellian, and the relative global entropy H g , the entropy relative to the global Maxwellian distribution, by Our goal here is to investigate numerically the long-time behavior of the solution f . If f is any reasonable solution of the Boltzmann equation, satisfying certain a priori bounds of compactness (in particular, ensuring that no kinetic energy is allowed to leak at large velocities), then it is possible to prove that f does indeed converge to the global Maxwellian distribution M g as t goes to +∞. Of course, obtaining these a priori bounds is extremely difficult; as a matter of fact, they have been established only in the spatially homogeneous situation (which means that the distribution function does not depend on the position variable x, see the survey in [49]) or in a close-to equilibrium setting (see in particular [25] for the torus), and it still constitutes a famous open problem for spatially inhomogeneous initial data far from equilibrium. More recently, Desvillettes and Villani [16], Guo and Strain [47] were interested in the study of rates of convergence for the full Boltzmann equation. Roughly speaking in [16], the authors proved that if the solution to the Boltzmann equation is smooth enough and satisfies bounds from below of the form (although this bound can be shown to be a consequence of the regularity bounds, see [30]) then (with constructive bounds) which means that the solution converges almost exponentially fast to the global equilibrium (namely with polynomial rate O(t −r ) with r as large as wanted).
The solution f to the Boltzmann equation satisfies the formula of additivity of the entropy: the entropy can be decomposed into the sum of a purely hydrodynamic part, and (by contrast) of a purely kinetic part. In terms of H functional: one can write In fact, we can also show that Moreover, the Csiszàr-Kullback-Pinsker inequality asserts that (when the total mass of the solution is normalized to 1) In other words, controlling the speed of convergence of the entropy to its equilibrium value is enough to control the speed of convergence of the solution to equilibrium, in very strong sense.
Moreover in [16], Desvillettes and Villani conjectured that time oscillations should occur on the evolution of the relative local entropy. In fact their proof does not rule out the possibility that the entropy production undergoes important oscillations in time, and actually most of the technical work is caused by this possibility.

Description and interpretation of the results.
Here, we performed simulations on the full Boltzmann equation in a simplified geometry (one dimension of space, two dimensions of velocity, periodic boundary conditions, fixed Knudsen number) with the fast spectral method to observe the evolution of the entropy and to check numerically if such oscillations occur. Clearly this test is challenging for a numerical method due to the high accuracy required to capture such oscillating behavior.
Then, we consider an initial datum as a perturbation of the global equilibrium M g for some constants A 0 > 0 and k 0 = 2π/L. In Figures 8 and 9, we are indeed able to observe oscillations in the entropy production and in the hydrodynamic entropy. The strength of the oscillations depends a lot on the length L of the domain, which is consistent with the fact that such oscillations are never observed in the spatially homogeneous case (L = 0). The superimposed curves yield the time evolution respectively of the total H functional and of its kinetic part. In all cases, a local Maxwellian distribution is chosen for initial datum; the first plot corresponds to L = 1 and the second one to L = 4. Some slight oscillations can be seen in the case L = 1, but what is most striking is that after a short while, the kinetic entropy is very close to the total entropy: an indication that the solution evolves basically in a spatially homogeneous way (contrary to the intuition of the hydrodynamic regime). On the contrary, in the case L = 4, the oscillations are much more important in frequency and amplitude (note that this is a logarithmic plot): the solution "hesitates" between states where it is very close to hydrodynamic, and states where it is not at all. Further note that the equilibration is much more rapid when the box is small, and that the convergence seems to be exponential.
It is in fact possible to give a simple interpretation of these oscillations thanks to the work [17]. Since this effect is observed near the global equilibrium one can replace the Boltzmann collision operator by the linearized Boltzmann collision operator (moreover the oscillation effect is effectively observed for the linearized Boltzmann collision operator as well). Then it is straightforward that the computations above correspond to observing the time evolution of one Fourier mode in x (here with frequency k 0 and amplitude A 0 ). Hence by an obvious rescaling, this evolution is given by the semi-groupT k 0 /L (t), whereT k (t) is defined in [17] (this is the semi-group for the k-th Fourier mode in x for the linearized equation). An asymptotic study of the spectrum of its infinitesimal generator for small frequencies k was done in [17]. The dominant term in terms of long time behavior (i.e., the one with the lower rate of decrease) is given by the (d + 2) "hydrodynamical eigenvalues". Moreover explicit computations are available for the expansions of these eigenvalues according to ε = |k| near k = 0.
At first order in ε, the eigenvalues vanish, except for two of them, which are purely imaginary. They are given by Therefore for |k 0 |/L << 1 (realized for instance when k 0 is fixed and L is large enough), this analysis gives us the dominant imaginary term in the eigenvalues. In this regime, one should thus observe oscillations with frequency 1 + 2/d |k 0 |/L. Thus the period of oscillations should be given by 2π(1 + 2/d) −1/2 L/|k 0 |, which can be checked with the numerical simulations. Indeed, in Table 4, we give the ratio of the period of oscillations ω with the length box L. The numerical results agree well with the analytical computations ω/L ≃ 1/ √ 2. We also observe that the damping rate is related to the length box and is proportional to 1/L 2 when L becomes large (see Table 4, α L 2 ≃ constant). This is coherent with the fact that no real value occurs in the "hydrodynamical" eigenvalues until the second order in ε = |k|. The coefficients for the order 2 in the expansion are computed in [17]; they are purely real and they can be expressed simply in terms of the dimension d, the viscosity coefficient η and the heat conductivity λ of the gas (indeed these coefficients are related with the Navier-Stokes limit of the Boltzmann equation). Namely they are given by Therefore for |k 0 |/L << 1, the damping rate is given by the minimum among these values.
To conclude these tests, we performed a last numerical experiment to evaluate the robustness of the theory of the trend to equilibrium. We have chosen an initial datum which is far from the equilibrium (6.4) with v 0 = (1/2, 1/2) and v th = √ 3/2. We present the time evolution of the relative entropies H l and H g in log scale and observe that initially the entropy is strongly decreasing and when the distribution function becomes close to a local equilibrium, some oscillations appear with the good frequency ω/L = 1/ √ 2 and damping rate α (see Figure 9).

Remarks:
1. Now numerical methods for the Boltzmann equation -such as the one presented in this paper-become able to provide very accurate simulations of the transient regime towards equilibrium with reasonable cost, even in the inhomogeneous case. This could be used to explore numerically the spectrum of the linearized Boltzmann collision operator (in the homogeneous case), or, more interestingly, the spectrum of the linearized Boltzmann collision operator together with the transport term in the inhomogeneous case. For instance the exponential rate of convergence is directly readable on the figures above, and provides a numerical estimation of the spectral gap (that is the real part of the first non-zero eigenvalue) of this operator (it is known since Ukai [48] that this operator has a spectral gap in the torus, see also [13]). Moreover by a frequency analysis of the curve of the time evolution of the relative entropy or the L 1 distance to the equilibrium, it could be possible also to describe other eigenvalues: as long as they have different imaginary part, it should be possible (in principle) to extract from the frequency analysis the curve corresponding to their contribution in the evolution semigroup, and then to compute their real part which corresponds to the exponential rate of decay of this curve.
Length box oscillation frequency ω ω/L damping rate α −α L 2 L = π/2 01.10 0.701 -6.521 16 Table 4. Influence of the length box: damping rate and oscillation frequency for the relative entropy with respect to the local Maxwellian H l (t) using 64 × 64 × 64 with A 0 = 0.1 in (6.3). [27] gave a detailed pointwise study of the Green function for the linearized Boltzmann equation in the domain x ∈ Ω = R. In particular in this case, there study shows that the long-time behavior is governed by "fluidlike waves" (corresponding to the waves of the linearized Euler and Navier-Stokes equations) whose amplitude decreases polynomially, whereas the amplitude of the "kinetic part" of the Green function decreases exponentially. We think it likely that this study could be extended to the torus, where the amplitude of the fluid and kinetic parts of the Green function should both decrease exponentially. Moreover the rate of decay of the kinetic part should not depend on the size of the box, whereas the rate of decay of the fluid part should do. Hence for a box small enough, the long-time behavior should be governed by the kinetic part of the Green function (that is like the spatially homogeneous Boltzmann equation), whereas for a box big enough, the long-time behavior should be governed by the "fluid-like waves". This is precisely what we observe numerically, and thus this theoretical study could provide a rigorous proof of the numerical observations above, at least in the linearized regime.

Conclusions
In this paper we have introduced and deeply tested a class of new fast algorithms for the computation of the Boltzmann collision operator. These methods allow to reduce the computational cost from O(n 2 ) to O(n log 2 n). We give computational evidence of the great performance of the schemes which can provide a dramatic speed up in computing time of deterministic schemes by making them competitive with Monte-Carlo methods where higher accuracy is required. A first numerical application to a non trivial problem in the space non homogeneous case confirms the strong computing potential of the new schemes.
Other methods such as singular value decomposition, fast multipole methods [9], separated representations in high-dimensional problems (see works by G.Beylkin for  instance [1]), wavelets, . . . could have been used to search for some decomposition of the form (4.19). However to our knowledge it is not known at now how to obtain the properties described above on this decomposition with these methods. H_l (t) H_g(t) y = -0.16 x -5.81 Figure 9. Influence of the length box: relative entropy with respect to the local Maxwellian H l (t) using 64 × 64 × 64 for L= 2π; 3π and 4π with A 0 = 0.1 in (6.4).
Villani for suggesting the numerical study of possible oscillations in the relaxation to equilibrium.