A new stability and convergence proof
of the Fourier-Galerkin spectral method
for the spatially homogeneous Boltzmann equation

The Boltzmann equation is one of the fundamental equations in kinetic theory and serves as a basic building block to connect microscopic Newtonian mechanics and macroscopic continuum mechanics [4, 22]. Albeit its wide applicability, numerical approximation of the Boltzmann equation is a challenging scientific problem due to the complicated structure of the equation (high-dimensional, nonlinear, and nonlocal). As such, the particle based direct simulation Monte Carlo method (DSMC) [2] has been widely used in various applications for its simplicity and low computational cost. Nevertheless, the stochastic method suffers from slow convergence and becomes extremely expensive when simulating non-steady and low-speed flows.

Since the pioneering work [18, 19], it has been realized that the Fourier-Galerkin spectral method offers a suitable framework to approximate the Boltzmann collision operator. First of all, it is a deterministic method and provides very accurate results compared with stochastic method. Secondly, the Boltzmann collision operator is translation-invariant and the Fourier basis exactly leverages this structure. Thirdly, after the Galerkin projection, the collision operator presents a convolution-like structure, which opens the possibility to further accelerate the method by the fast Fourier transform (FFT) [16, 9]. Because of the above reasons, over the past decade, the Fourier spectral method has become a very popular deterministic method for solving the Boltzmann equation and related collisional kinetic models, see for instance, [21, 8, 7, 14, 15], or the recent review article [5].

As opposed to its practical success, the theoretical study of the Fourier spectral method is quite limited, largely because the spectral approximation destroys the positivity of the solution, yet the positivity is one of the key properties to study the well-posedness of the equation. In [20], a positivity-preserving filter is applied to the equation to enforce the positivity of the solution. As a result, the stability of the method can be easily proved. However, the filter often comes with the price of significantly smearing the solution (hence destroying the spectral accuracy) and should be used only when the solution contains discontinuities (to suppress the oscillations caused by Gibbs phenomenon). Recently, a stability proof for the original Fourier spectral method is established in [6], where the authors provide a quite complete study of the method including both finite and long time behavior. The key strategy in [6] is to use the “spreading” or “mixing” property of the collision operator to show that the solution will become everywhere positive after a small time. Motivated by this work, we present in this paper a different well-posedness and stability proof. The main difference from [6] lies in that, instead of requiring the solution to be positive everywhere which is a stronger condition to achieve, we show that the $L^{2}$ norm of the negative part of the solution can be controlled as long as it is small initially. In other words, the solution is allowed to be negative for the method to remain stable. Therefore, our strategy does not rely on any sophisticated property of the collision operator and provides a simpler proof. In addition, we quantify clearly the requirement on the initial condition for the method to be stable, which includes both continuous and discontinuous functions.

We mention another line of research which develops the conservative-spectral approximation for the Boltzmann equation [10]. Apart from apparent differences (the Fourier-Galerkin method considered in this paper is based on domain truncation and periodization, while the method [10] is based on Fourier transform and no periodization is performed), a conservation subroutine is added to restore the mass, momentum, and energy conservation. As a consequence, the method is able to preserve the Maxwellian distribution as time goes to infinity. The stability and convergence of the method is recently established in [1], where the Fourier projection is only applied to the gain part of the collision operator. In contrast, both gain and loss terms are projected in our method, hence the loss term does not possess a definite sign.

The paper is essentially self-contained. In Section 2, we briefly review the Fourier-Galerkin spectral method for the spatially homogeneous Boltzmann equation. After that, we discuss the basic assumptions (e.g., the collision kernel and truncation parameters) used throughout the paper. The assumptions on the initial condition are addressed in Section 2.1, which will play an important role in proving the main result. In Section 3 (and Appendix), we provide some preliminary estimates on the truncated collision operator. These are known results in the whole space but some subtle differences appear in the torus. Section 4 presents our main result. We first conduct a $L^{2}$ estimate of the negative part of the solution and then prove a local existence/uniqueness result. Finally, the well-posedness and stability of the method on an arbitrary bounded time interval is established in Section 4.3 (Theorem 4.4). Facilitated with the stability result, the paper is concluded in Section 5 with a straightforward convergence and spectral accuracy proof of the method.

In this section, we review the Fourier-Galerkin spectral method for the spatially homogeneous Boltzmann equation. The presentation follows the formulation originally proposed in [19] which is the basis for many fast algorithms developed recently [9, 12, 13]. Here we limit the description to the extent that is sufficient for the following proof. At the end of the section, we discuss the basic assumptions used throughout the rest of the paper, in particular, the assumptions on the initial condition.

The spatially homogeneous Boltzmann equation reads

$$\partial_{t}f=Q(f,f),\quad t>0,\ v\in\mathbb{R}^{d},\ d\geq 2,$$

(2.1)

where $f=f(t,v)$ is the probability density function of time $t$ and velocity $v$, $Q$ is the collision operator describing the binary collisions among particles, whose bilinear form is given by

$$Q(g,f)(v)=\int_{\mathbb{R}^{d}}\int_{\mathbb{S}^{d-1}}B(|v-v_{*}|,\cos\theta)[g(v_{*}^{\prime})f(v^{\prime})-g(v_{*})f(v)]\,\mathrm{d}{\sigma}\,\mathrm{d}{v_{*}}.$$

(2.2)

In (2.2), $\sigma$ is a vector varying over the unit sphere $\mathbb{S}^{d-1}$, $v^{\prime}$ and $v_{*}^{\prime}$ are defined as

$$v^{\prime}=\frac{v+v_{*}}{2}+\frac{|v-v_{*}|}{2}\sigma,\quad v_{*}^{\prime}=\frac{v+v_{*}}{2}-\frac{|v-v_{*}|}{2}\sigma,$$

(2.3)

and $B\geq 0$ is the collision kernel. In this paper we will consider the kernel of the form

$$B(|v-v_{*}|,\cos\theta)=\Phi(|v-v_{*}|)b(\cos\theta),\quad\cos\theta=\frac{\sigma\cdot(v-v_{*})}{|v-v_{*}|},$$

(2.4)

whose kinetic part $\Phi$ is a non-negative function and angular part $b$ satisfies the Grad’s cut-off assumption

$$\int_{\mathbb{S}^{d-1}}b(\cos\theta)\,\mathrm{d}{\sigma}<\infty.$$

(2.5)

To apply the Fourier-Galerkin spectral method, we consider an approximated problem of (2.1) on a torus $\mathcal{D}_{L}=[-L,L]^{d}$:

$$\left\{\begin{split}&\partial_{t}f=Q^{R}(f,f),\quad t>0,\ v\in\mathcal{D}_{L},\\ &f(0,v)=f^{0}(v),\end{split}\right.$$

(2.6)

where the initial condition $f^{0}$ is a non-negative periodic function, $Q^{R}$ is the truncated collision operator defined by

$$\begin{split}Q^{R}(g,f)(v)&=\int_{\mathcal{B}_{R}}\int_{\mathbb{S}^{d-1}}\Phi(|q|)b(\sigma\cdot\hat{q})\left[g(v_{*}^{\prime})f(v^{\prime})-g(v-q)f(v)\right]\,\mathrm{d}{\sigma}\,\mathrm{d}{q}\\ &=\int_{\mathbb{R}^{d}}\int_{\mathbb{S}^{d-1}}\mathbf{1}_{|q|\leq R}\Phi(|q|)b(\sigma\cdot\hat{q})\left[g(v_{*}^{\prime})f(v^{\prime})-g(v-q)f(v)\right]\,\mathrm{d}{\sigma}\,\mathrm{d}{q},\end{split}$$

(2.7)

where a change of variable $v_{*}\rightarrow q=v-v_{*}$ is applied and the new variable $q$ is truncated to a ball $\mathcal{B}_{R}$ with radius $R$ centered at the origin. We write $q=|q|\hat{q}$ with $|q|$ being the magnitude and $\hat{q}$ being the direction. Accordingly,

$$v^{\prime}=v-\frac{q-|q|\sigma}{2},\quad v_{*}^{\prime}=v-\frac{q+|q|\sigma}{2}.$$

(2.8)

In practice, the values of $L$ and $R$ are often chosen by an anti-aliasing argument [19]: assume that $\text{Supp}(f^{0}(v))\subset\mathcal{B}_{S}$, then one can take

$$R=2S,\quad L\geq\frac{3+\sqrt{2}}{2}S.$$

(2.9)

Given an integer $N\geq 0$, we then seek a truncated Fourier series expansion of $f$ as

$$f(t,v)\approx f_{N}(t,v)=\sum\limits_{k=-N/2}^{N/2}f_{k}(t)\mathrm{e}^{\mathrm{i}\frac{\pi}{L}k\cdot v}\in\mathbb{P}_{N},$$

(2.10)

where

$$\mathbb{P}_{N}=\text{span}\left\{\mathrm{e}^{\mathrm{i}\frac{\pi}{L}k\cdot v}\Big{|}-N/2\leq k\leq N/2\right\},$$

(2.11)

equipped with inner product

$$\langle f,g\rangle=\frac{1}{(2L)^{d}}\int_{\mathcal{D}_{L}}f\bar{g}\,\mathrm{d}v.$$

(2.12)

Substituting $f_{N}$ into (2.6) and conducting the Galerkin projection onto the space $\mathbb{P}_{N}$ yields

$$\left\{\begin{split}&\partial_{t}f_{N}=\mathcal{P}_{N}Q^{R}(f_{N},f_{N}),\quad t>0,\ v\in\mathcal{D}_{L},\\ &f_{N}(0,v)=f_{N}^{0}(v),\end{split}\right.$$

(2.13)

where $\mathcal{P}_{N}$ is the projection operator: for any function $g$,

$$\mathcal{P}_{N}g=\sum_{k=-N/2}^{N/2}\hat{g}_{k}\mathrm{e}^{\mathrm{i}\frac{\pi}{L}k\cdot v},\quad\hat{g}_{k}=\langle g,\mathrm{e}^{\mathrm{i}\frac{\pi}{L}k\cdot v}\rangle,$$

(2.14)

$f_{N}^{0}\in\mathbb{P}_{N}$ is the initial condition to the numerical system and should be a reasonable approximation to $f^{0}$. More discussion on the initial condition will be given in Section 2.1, which in fact plays an important role in the following proof.

Writing out each Fourier mode of (2.13), we obtain

$$\left\{\begin{split}&\partial_{t}f_{k}=Q^{R}_{k},\quad-N/2\leq k\leq N/2,\\ &f_{k}(0)=f^{0}_{k},\end{split}\right.$$

(2.15)

with

$$Q_{k}^{R}:=\langle Q^{R}(f_{N},f_{N}),\mathrm{e}^{\mathrm{i}\frac{\pi}{L}k\cdot v}\rangle,\quad f^{0}_{k}:=\langle f_{N}^{0},\mathrm{e}^{\mathrm{i}\frac{\pi}{L}k\cdot v}\rangle.$$

(2.16)

Using the definition in (2.7) and orthogonality of the Fourier basis, we can derive that

$$Q_{k}^{R}=\sum\limits_{\begin{subarray}{c}l,m=-N/2\\ l+m=k\end{subarray}}^{N/2}G(l,m)f_{l}f_{m},$$

(2.17)

where the weight $G$ is given by

$$\begin{split}G(l,m)&=\int_{\mathcal{B}_{R}}\int_{\mathbb{S}^{d-1}}\Phi(|q|)b(\sigma\cdot\hat{q})\left[\mathrm{e}^{-\mathrm{i}\frac{\pi}{2L}(l+m)\cdot q+\mathrm{i}\frac{\pi}{2L}|q|(l-m)\cdot\sigma}-\mathrm{e}^{-\mathrm{i}\frac{\pi}{L}m\cdot q}\right]\,\mathrm{d}\sigma\,\mathrm{d}q\\ &=\int_{\mathcal{B}_{R}}\mathrm{e}^{-\mathrm{i}\frac{\pi}{L}m\cdot q}\left[\int_{\mathbb{S}^{d-1}}\Phi(|q|)b(\sigma\cdot\hat{q})(\mathrm{e}^{\mathrm{i}\frac{\pi}{2L}(l+m)\cdot(q-|q|\sigma)}-1)\,\mathrm{d}\sigma\right]\,\mathrm{d}q.\end{split}$$

(2.18)

The second equality above is obtained by switching two variables $\sigma\leftrightarrow\hat{q}$ in the gain part of $G(l,m)$. In the direct Fourier spectral method, $G(l,m)$ is precomputed since it is independent of the solution. Then in the online computation, the sum (2.17) is evaluated directly.

Note that the solution $f$ to the original problem (2.6) is always non-negative which is the key to many stability estimates. However, the solution $f_{N}$ to the numerical system (2.13) is not necessarily non-negative due to the spectral projection which constitutes the main difficulty in the numerical analysis. Luckily, by virtue of the Fourier spectral method, mass is always conserved which provides some control of the solution. Precisely, we have

We now introduce some assumptions and notations that will be used throughout the rest of this paper.

• (1)

  The truncation parameters $L$ and $R$ in (2.6) satisfy

  $$L\geq R>0.$$

  (2.22)

  Note that the choice (2.9) implies $L\geq(3+\sqrt{2})R/4$ hence the above condition is satisfied.

• (2)

  The kinetic part of the collision kernel (2.4) satisfies

  $$\left\|\mathbf{1}_{|v|\leq R}\Phi(|v|)\right\|_{L^{\infty}(\mathcal{D}_{L})}<\infty.$$

  (2.23)

  Note that all power law hard potentials $\Phi(|v|)=|v|^{\gamma}$ ($0\leq\gamma\leq 1$) as well as the “modified” soft potentials $\Phi(|v|)=(1+|v|)^{\gamma}$ ($-d<\gamma<0$) satisfy this condition.

• (3)

  The angular part of the collision kernel (2.4) has been replaced by its symmetrized version⁵⁵5This symmetrization can readily reduce the computational cost by a half (integration over the whole sphere is reduced to half sphere) so it also has important implications for numerical purpose, see [9].:

  $$\left[b(\cos\theta)+b(\cos\left(\pi-\theta\right))\right]\mathbf{1}_{0\leq\theta\leq\pi/2},$$

  (2.24)

  and satisfies the cut-off assumption (2.5).

For a periodic function $f(v)$ in $\mathcal{D}_{L}$, we define its Lebesgue norm and Sobolev norm as follows:

$$\|f\|_{L^{p}_{\text{per}}(\mathcal{D}_{L})}=\left(\int_{\mathcal{D}_{L}}|f(v)|^{p}\,\mathrm{d}{v}\right)^{1/p},\quad\|f\|_{{H^{k}_{\text{per}}}(\mathcal{D}_{L})}=\left(\sum_{|\nu|\leq k}\|\partial_{v}^{\nu}f\|_{L^{2}_{\text{per}(\mathcal{D}_{L})}}^{2}\right)^{1/2},$$

(2.25)

where $k\geq 0$ is an integer and $\nu$ is a multi-index. “per” indicates the function is periodic and will not be included in the following for simplicity.

Except in Section 3, we do not track explicitly the dependence of constants on the truncation parameters $R$, $L$, dimension $d$, and the collision kernel $B$.

For a function $f(v)$ in $\mathcal{D}_{L}$, we define its positive and negative parts as

$$f^{+}(v)=\max\limits_{v\in\mathcal{D}_{L}}\{f(v),0\},\quad f^{-}(v)=\max\limits_{v\in\mathcal{D}_{L}}\{-f(v),0\},$$

(2.26)

so that $f=f^{+}-f^{-}$ and $|f|=f^{+}+f^{-}$.

In this section, we prove some important estimates for the truncated collision operator (2.7). Since its gain term and loss term possess quite different properties, we consider

$$\begin{split}Q^{R,+}(g,f)(v)&:=\int_{\mathbb{R}^{d}}\int_{\mathbb{S}^{d-1}}\mathbf{1}_{|q|\leq R}\Phi(|q|)b(\sigma\cdot\hat{q})g(v_{*}^{\prime})f(v^{\prime})\,\mathrm{d}{\sigma}\,\mathrm{d}{q},\\ Q^{R,-}(g,f)(v)&:=\int_{\mathbb{R}^{d}}\int_{\mathbb{S}^{d-1}}\mathbf{1}_{|q|\leq R}\Phi(|q|)b(\sigma\cdot\hat{q})g(v-q)f(v)\,\mathrm{d}{\sigma}\,\mathrm{d}{q},\end{split}$$

(3.1)

separately whenever appropriate.

In this section, we establish the well-posedness and stability of the Fourier-Galerkin spectral method (2.13) on an arbitrary bounded time interval $[0,T]$. The main strategy of the proof is as follows: In Section 4.1 we prove some $L^{2}$ and $H^{k}$ estimates of the solution under a priori $L^{1}$ bound of $f_{N}$, among which the key result is the $L^{2}$ estimate of the negative part of the solution (Proposition 4.2). Proposition 4.3 is a local existence and uniqueness result over a small time interval $[t_{0},t_{0}+\tau]$. Finally the main result is presented in Theorem 4.4, where we show that when $N$ is large enough the negative part of the solution can be controlled over time $[0,\tau]$. Due to mass conservation, this consequently implies the initial $L^{1}$ bound of the solution can be restored at time $\tau$. Therefore, we can repeat the procedure iteratively to build the solution up to final time $T$ (the estimates on $N$ and $\tau$ are done carefully at the beginning so that the same values can be used in the following iteration).