Approximation to Singular Quadratic Collision Model in Fokker-Planck-Landau Equation

TL;DR

This paper introduces a Hermite-Galerkin spectral method for efficiently solving the Fokker-Planck-Landau equation with a singular quadratic collision model, accurately capturing moments despite singularities.

Contribution

It presents a novel spectral approach combining Burnett polynomials and Hermite expansion to handle singular collision models in kinetic equations.

Findings

Accurately captures low-order moments.

Handles very singular collision models smoothly.

Demonstrates satisfactory accuracy and performance.

Abstract

We propose a Hermite-Galerkin spectral method to numerically solve the spatially homogeneous Fokker-Planck-Landau equation with singular quadratic collision model. To compute the collision model, we adopt a novel approximation formulated by a combination of a simple linear term and a quadratic term very expensive to evaluate. Using the Hermite expansion, the quadratic term is evaluated exactly by calculating the spectral coefficients. To deal with singularities, we make use of Burnett polynomials so that even very singular collision model can be handled smoothly. Numerical examples demonstrate that our method can capture low-order moments with satisfactory accuracy and performance.

Equations

$$\frac{\partial f}{\partial t}+{\nabla }_{\mathbf{x}}\cdot (\mathbf{v}f)=\mathcal{Q}[f],t\in {\mathbb{R}}^{+},\mathbf{x}\in {\mathbb{R}}^3,\mathbf{v}\in {\mathbb{R}}^3,$$

$$\mathcal{Q}[f](t,\mathbf{x},\mathbf{v})={\nabla }_{\mathbf{v}}\cdot {\int }_{{\mathbb{R}}^3}A(\mathbf{v}-{\mathbf{v}}_{\ast })(f({\mathbf{v}}_{\ast }){\nabla }_{\mathbf{v}}f(\mathbf{v})-f(\mathbf{v}){\nabla }_{{\mathbf{v}}_{\ast }}f({\mathbf{v}}_{\ast }))\mathrm{d}{\mathbf{v}}_{\ast },$$

$$A(\mathbf{v})=\Psi (|\mathbf{v}|)\Pi (\mathbf{v}),$$

$$\Psi (\mathbf{v}):=\Lambda |\mathbf{v}{|}^{\gamma +2},$$

$$\frac{\partial f}{\partial t}=\mathcal{Q}[f],t\in {\mathbb{R}}^{+},\mathbf{v}\in {\mathbb{R}}^3.$$

$${\mathcal{M}}_{\rho ,\mathbf{u},\theta }(\mathbf{v}):=\frac{\rho }{(2\pi \theta {)}^{3/2}}\exp \left(-\frac{|\mathbf{v}-\mathbf{u}{|}^2}{2\theta }\right),$$

$$\rho ={\int }_{{\mathbb{R}}^3}f(t,\mathbf{v})\mathrm{d}\mathbf{v},\mathbf{u}=\frac{1}{\rho }{\int }_{{\mathbb{R}}^3}\mathbf{v}f(t,\mathbf{v})\mathrm{d}\mathbf{v},\theta =\frac{1}{3\rho }{\int }_{{\mathbb{R}}^3}|\mathbf{v}-\mathbf{u}{|}^{2}f(t,\mathbf{v})\mathrm{d}\mathbf{v}.$$

$$q_i=\frac{1}{2}{\int }_{{\mathbb{R}}^3}|\mathbf{v}-\mathbf{u}{|}^2(v_i-u_i)f\mathrm{d}\mathbf{v},i=1,2,3,$$

$${\sigma }_{ij}={\int }_{{\mathbb{R}}^3}\left((v_i-u_i)(v_j-u_j)-\frac{1}{3}{\delta }_{ij}|\mathbf{v}-\mathbf{u}{|}^2\right)f\mathrm{d}\mathbf{v},i,j=1,2,3.$$

$$\rho =1,\mathbf{u}=0,\theta =1,$$

$$\mathcal{M}(\mathbf{v}):=\frac{1}{(2\pi {)}^{3/2}}\exp \left(-\frac{|\mathbf{v}{|}^2}{2}\right).$$

$$q_i=\frac{1}{2}{\int }_{{\mathbb{R}}^3}|\mathbf{v}{|}^2{v}_{i}f\mathrm{d}\mathbf{v},{\sigma }_{ij}={\int }_{{\mathbb{R}}^3}\left(v_i{v}_j-\frac{1}{3}{\delta }_{ij}|\mathbf{v}{|}^2\right)f\mathrm{d}\mathbf{v},i,j=1,2,3.$$

$${\mathcal{Q}}^{\mathrm{linear}}[f]=(D-1){\nabla }_{\mathbf{v}}\cdot (\nabla f+f\mathbf{v}),$$

$$f(t,\mathbf{v})=\sum _{|\alpha |=0}^{+\infty }{f}_{\alpha }(t)H^{\alpha }(\mathbf{v})\mathcal{M}(\mathbf{v}),$$

$$H^{\alpha }(\mathbf{v})=\frac{(-1{)}^n}{\mathcal{M}(\mathbf{v})}\frac{{\partial }^{|\alpha |}}{\partial v_1^{{\alpha }_1}\partial v_2^{{\alpha }_2}\partial v_3^{{\alpha }_3}}\mathcal{M}(\mathbf{v}),$$

$${\int }_{{\mathbb{R}}^3}{H}^{\alpha }(\mathbf{v})H^{\beta }(\mathbf{v})\mathcal{M}(\mathbf{v})\mathrm{d}\mathbf{v}={\delta }_{\alpha ,\beta }\alpha !,$$

$$f_{e_i}=0,i=1,2,3,\sum _{i=1}^3{f}_{2e_i}=0,$$

$$q_j=2f_{3e_j}+\sum _{k=1}^3{f}_{e_j+2e_k},{\sigma }_{ij}=(1+{\delta }_{ij})f_{e_i+e_j},i,j=1,2,3.$$

$${\mathcal{Q}}^{\mathrm{linear}}[f]=\sum _{|\alpha |=0}^{+\infty }{Q}_{\alpha }^{\mathrm{linear}}{H}^{\alpha }(\mathbf{v})\mathcal{M}(\mathbf{v}),Q_{\alpha }^{\mathrm{linear}}=-(D-1)|\alpha |f_{\alpha }.$$

$$\mathcal{Q}[f](\mathbf{v})=\sum _{|\alpha |=0}^{+\infty }{Q}_{\alpha }{H}^{\alpha }(\mathbf{v})\mathcal{M}(\mathbf{v}).$$

$$Q_{\alpha }=\frac{1}{\alpha !}\int H^{\alpha }(\mathbf{v})\mathcal{Q}[f](\mathbf{v})\mathrm{d}\mathbf{v}=\sum _{|\lambda |=0}^{+\infty }\sum _{|\kappa |=0}^{+\infty }{A}_{\alpha }^{\lambda ,\kappa }{f}_{\lambda }{f}_{\kappa },$$

$$A_{\alpha }^{\lambda ,\kappa }=\frac{1}{\alpha !}{\int }_{{\mathbb{R}}^3}{H}^{\alpha }(\mathbf{v}){\nabla }_{\mathbf{v}}\cdot {\int }_{{\mathbb{R}}^3}A(\mathbf{v}-{\mathbf{v}}_{\ast })$$

$$\displaystyle\qquad\Big{(}H^{{\lambda}}(\boldsymbol{v}_{\ast})\mathcal{M}(\boldsymbol{v}_{\ast})\nabla_{\boldsymbol{v}}\left(H^{{\kappa}}(\boldsymbol{v})\mathcal{M}(\boldsymbol{v})\right)-H^{{\lambda}}(\boldsymbol{v})\mathcal{M}(\boldsymbol{v})\nabla_{\boldsymbol{v}_{\ast}}\left(H^{{\kappa}}(\boldsymbol{v}_{\ast})\mathcal{M}(\boldsymbol{v}_{\ast})\right)\Big{)}\,\mathrm{d}\boldsymbol{v}_{\ast}\,\mathrm{d}\boldsymbol{v}.$$

$$A_{\alpha }^{\lambda ,\kappa }=$$

$$q^{[s]}=\alpha -e_s-p,r^{[t]}=\lambda +\kappa +e_t-p,a_{p,q}^{\lambda ,\kappa }=\prod _{i=1}^3{a}_{p_i{q}_i}^{{\lambda }_i{\kappa }_i},s,t=1,2,3.$$

$$a_{pq}^{\lambda \kappa }=2^{-(p+q)/2}\lambda !\kappa !\sum _{s=\max (0,p-\kappa )}^{\min (p,\lambda )}\frac{(-1{)}^{q-\lambda +s}}{s!(\lambda -s)!(p-s)!(q-\lambda +s)!},$$

$$B_p^q(\gamma ,s,t):=-G_{st}(\gamma ,p,q)+{\delta }_{st}\sum _{r=1}^3{G}_{rr}(\gamma ,p,q),$$

$$G_{st}(\gamma ,p,q)={\int }_{\mathbf{g}\in {\mathbb{R}}^3}|\mathbf{g}{|}^{\gamma }{g}_s{g}_t{H}^p(\mathbf{g})H^q(\mathbf{g})\mathcal{M}(\mathbf{g})\mathrm{d}\mathbf{g},s,t=1,2,3.$$

$$B_{\hat{\alpha }}(\mathbf{v})=\sqrt{\frac{2^{1-{\hat{\alpha }}_1}{\pi }^{3/2}{\hat{\alpha }}_3!}{\Gamma ({\hat{\alpha }}_3+{\hat{\alpha }}_1+3/2)}}{L}_{{\hat{\alpha }}_3}^{({\hat{\alpha }}_1+1/2)}\left(\frac{|\mathbf{v}{|}^2}{2}\right)|\mathbf{v}{|}^{{\hat{\alpha }}_1}{Y}_{{\hat{\alpha }}_1}^{{\hat{\alpha }}_2}\left(\frac{\mathbf{v}}{|\mathbf{v}|}\right),$$

$$\hat{\alpha }=({\hat{\alpha }}_1,{\hat{\alpha }}_2,{\hat{\alpha }}_3{)}^T,{\hat{\alpha }}_1,{\hat{\alpha }}_3\in \mathbb{N},{\hat{\alpha }}_2=-{\hat{\alpha }}_1,\cdots ,{\hat{\alpha }}_1.$$

$$L_n^{(\beta )}(x)=\frac{x^{-\beta }\exp (x)}{n!}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n+\beta }\exp (-x)\right],$$

Full text

Keywords: Fokker-Planck-Landau equation; Hermite-Galerkin spectral method; Burnett polynomials; Quadratic collision operator; Super singularity

1 Introduction

The Fokker-Planck-Landau (FPL) equation is a common kinetic model in plasma physics, accelerator physics and astrophysics. It describes binary collision between charged particles with long-range Coulomb interaction, and is represented by a nonlinear partial integro-differential equation.

As a classical result, the FPL operator is the limit of the Boltzmann operator for a sequence of scatting cross sections which converge in a convenient sense to a delta function at zero scattering angle [8]. The original derivation of the equation based on this idea is due to Landau [22], and then several work has been devoted to this problem, such as [2, 12, 29]. Recently Villani [33] has obtained a rigorous proof of this asymptotic problem in the space homogeneous scenario. For the mathematical properties of the FPL equation, such as the existence of the solutions, we refer the reader to Villani [32] and the reference therein.

The numerical solution of the nonlinear kinetic equations, such as FPL equation, also represents a real challenge for numerical method. This is essentially due to the non-linearity, as well as the high dimension of variables, which is seven for the full problem. Moreover, the complex three dimensional integro-differential stiff advection-diffusion operator in velocity space is also remarkably difficult to deal with due to the high singularity. Besides, this integration has to be handled carefully since it is closely related to the macroscopic properties, for example the collision term does not change the total mass, momentum and energy. Several numerical approaches have been brought up to solve FPL equation. Generally speaking, there are two kinds of methods, stochastic methods and the deterministic methods. For the stochastic methods, DSMC method which is widely in the simulation of Boltzmann equation[4] is adopted to solve FPL equation. A detail discussion about the stochastic method is beyond the scope of this paper, and we refer the reader to [15, 5] for a much more complete treatment. For the deterministic methods, due to the complex form of the FPL operator, several numerical approaches are devoted to the simpler diffusive Fokker-Planck model [17, 37], the space homogeneous situations in the isotropic case [6] or cylindrically symmetric problems [26]. Moreover, Villani [34] has brought up a linear collision model for the Maxwell molecules. The construction of conservative and entropy schemes for the space homogeneous case has been proposed in [13, 7], where the main physical properties are all satisfied. But the direct implementation of such schemes are all quite expensive. Several fast approximated algorithms to reduce the complexity of these methods, based on multipole expansion [23] or multigrid techniques [6] have been proposed. A fast spectral method based on Fourier spectral approximation of the collision operator is introduced in [25], and it is then also utilized to solve inhomogeneous FPL equation [36, 16]. For the numerical stiffness of the Fokker-Planck collision operator, the implicit time scheme is also studied [24, 31]. There is a certain kind of asymptotic-preserving method that seeks to accelerate the solution of the FPL equation by the so-called penalization techniques [21].

As another kind of spectral method, Hermite spectral method is also utilized to solve FPL equation. Hermite method, where the basis functions with weighted orthogonality in $\mathbb{R}^{3}$ are employed, dates back to Grad’s work [19] where it is used to solve the Boltzmann equation and is known as the moment method ever since. Besides, the expansion with respect to Burnett polynomials was proposed in [9, 18] to find the coefficients of the collision term of the expansion in the Hermite basis. Using the Hermite expansion, it is still a tough job to evaluate the exact coefficients in the expansion of the collision operator, since the computational cost for the quadratic from is hardly bearable and novel models need to be introduced. In a recent work [35], the explicit expressions of all the coefficients in the Hermite spectral method for the quadratic Boltzmann collision operator are presented, and the new collision model which can preserve the physical properties and reduce computational cost at the same time was brought up using these coefficients. It is much harder to evaluate these coefficients for the quadratic FPL collision operator compared to the Boltzmann equation, because of the high singularity and the operator of partial derivative. In [27], the coefficients for the Coulombian case were evaluated numerically, and the explicit form was listed for the first few moments.

Inspired by these work, we in this paper are devoted to the numerical method for FPL equation with quadratic collision model, which may be very singular. Following the approach in [35], we approximate the collision model as the combination of a simple linear term and a quadratic term. The idea is to take only a portion in the truncated series expansion to be treated “quadratically”, and the remaining part is approximated by the linear collision operator brought up by Villani [34]. This may greatly reduce the computational cost and we can still capture the evolution of physical variables accurately. The linear term can be handled easily, while the difficulty imposed by the singularity in the quadratic collision model remains. We reveal that by making use of Burnett polynomials, the singular part of the integral in the collision operator can be handled smoothly. For the typical case that the repulsive force between molecules is proportional to a negative power of their distance, our method can handle problems where the index for the power of distance is as great as $-5$, in comparison to the index fixed as $-3$ in [27]. To deal with the remaining part in the quadratic term without singularity, the Hermite-Galerkin spectral method is then adopted. We derive the explicit formulae for all the coefficients in the Hermite expansion of the collision operator, and these formulae can all be evaluated exactly offline for immediate applications. Thus eventually the quadratic term is able to be evaluated efficiently.

The rest of this paper is organized as follows. In Section 2, we briefly review the FPL equation and the Hermite expansion of the distribution function. In Section 3, we first give an explicit expression of the series expansion of the quadratic collision operator and then introduce precisely how to deal with the singularity by Burnett polynomials. The construction of the approximated collision model is presented in Section 4. Some numerical experiments verifying the effectiveness of our methods are carried out in Section 5. The concluding remarks and detailed derivation of the expansions are given in Section 6 and 7 respectively.

2 FPL equation and Hermite expansion

We will first give a brief review of the FPL equation, and then introduce the Hermite spectral method for the expansion of the distribution function.

2.1 FPL equation

The Fokker-Planck-Landau equation is a prevalent kinetic model in plasma physics, describing the state of the particles in terms of a distribution function $f(t,\boldsymbol{x},\boldsymbol{v})$, where $t$ is the time coordinate, $\boldsymbol{x}$ represents the spatial coordinates, and $\boldsymbol{v}$ stands for the velocity of particles. The governing equation of $f$ is

[TABLE]

where $\mathcal{Q}[f]$ is the collision operator with a quadratic form

[TABLE]

where $A$ depends on the interaction between particles and is a $3\times 3$ negative and symmetric matrix in the form [25] of

[TABLE]

where $\Psi$ is a non-negative radial function, and $\Pi(\boldsymbol{v})$ is the orthogonal projection upon the space orthogonal to $\boldsymbol{v}$, as $\Pi_{ij}(\boldsymbol{v})=\delta_{ij}-\frac{v_{i}v_{j}}{|\boldsymbol{v}|^{2}}$.

We are primarily concerned with the IPL model, for which the force between two molecules is always repulsive and proportional to a negative power of their distance. In this case, the function $\Psi(\boldsymbol{v})$ has the form

[TABLE]

where $\Lambda>0$ is a constant and $\gamma$ is the index of the power of distance. This equation is obtained as a limit of the Boltzmann equation, when all the collisions become grazing [14]. In the case of the Boltzmann equation, different $\gamma$ lead to different models. The case $\gamma>0$ corresponds to the “hard potential” case, whereas for $\gamma<0$, it corresponds to the case of “soft potential”. In the critical case $\gamma=0$, the gas molecules are referred to as “Maxwell molecules”. Another case of interest is when $\gamma=-3$ of the Coulombian case, which is a very important model for applications in plasma.

We shall focus on the numerical approximation of $\mathcal{Q}[f]$, especially when $\gamma$ is very small. Our model of approximating the collision operator is best illustrated in the spatially homogeneous FPL equation case, namely

[TABLE]

As a classical result in kinetic equations, the steady state solution of this equation takes the form of the Maxwellian:

[TABLE]

where the density $\rho$, velocity $\boldsymbol{u}$ and temperature $\theta$ are defined as follows

[TABLE]

Moreover, the physical variables such as the heat flux $q_{i}$ and the stress tensor $\sigma_{ij}$ are also of interest. They are defined as

[TABLE]

Similar to the Boltzmann equation, the collision operator preserves in time the macroscopic quantities mass, momentum and energy. Therefore, those are invariant quantities under evolution, and (2.7) holds for any $t$. Thus we can obtain

[TABLE]

by selecting proper frame of reference and applying appropriate non-dimensionalization. Now the Maxwellian (2.6) is simply reduced to

[TABLE]

The heat flux and stress tensor are reduced into

[TABLE]

The normalization (2.8) shall always be assumed in the following context.

In the literature, the complicated form of the collision operator $\mathcal{Q}[f]$ is handled by introducing approximations of less complexity. For instance, for the Maxwell molecules with $\Lambda=1$, if the distribution function $f$ is radially symmetric, which is a property to be preserved under time evolution, the collision operator can be rewritten as

[TABLE]

which was proposed by C. Villani [34] . Here $D$ is the dimension of the velocity space, and we always set $D=3$ in the context. In this case, the FPL equation is reduced into the linear Fokker-Planck equation (FP), which can be used to describe the relation of Brownian molecules in a gas.

Due to the complex form of the FPL operator, several numerical approaches are devoted to the simpler diffusive Fokker-Planck model or on the reduced collision models [28, 3]. Hence it is of high necessity to develop efficient numerical methods for the original FPL equation with quadratic collision operator.

2.2 Series expansion of distribution function

Our numerical discretization shall be based on the series expansion in the weighted $L^{2}$ space of the distribution function $\mathcal{F}=L^{2}(\mathbb{R}^{3};\mathcal{M}^{-1}\,\mathrm{d}\boldsymbol{v})$:

[TABLE]

where $\mathcal{M}(\boldsymbol{v})$ is the Maxwellian, and ${\alpha}=(\alpha_{1},\alpha_{2},\alpha_{3})^{T}$ is a three-dimensional multi-index, and $|{\alpha}|=\alpha_{1}+\alpha_{2}+\alpha_{3}$. In (2.11), $H^{{\alpha}}(\boldsymbol{v})$ are the Hermite polynomials defined as follows:

Definition 1** (Hermite polynomials).**

For $\alpha_{i}\in\mathbb{N},i=1,2,3$, the Hermite polynomial $H^{{\alpha}}(\boldsymbol{v})$ is defined as

[TABLE]

where $\mathcal{M}(\boldsymbol{v})$ is given in (2.9).

The expansion (2.2) was introduced to solve Boltzmann equations [19], where such an expansion was invoked to deploy moment methods. We can derive moments based on the coefficients $f_{{\alpha}}$ from the orthogonality of Hermite polynomials

[TABLE]

where $\delta_{{\alpha},{\beta}}$ is defined as $\delta_{{\alpha},{\beta}}=\prod\limits_{i=1}^{3}\delta_{\alpha_{i},\beta_{i}}$ and ${\alpha}!=\prod\limits_{i=1}^{3}\alpha_{i}!$. For example, by the orthogonality aforementioned, we can insert the expansion (2.11) into the definition of $\rho$ in (2.7) to get $f_{\boldsymbol{0}}=\rho$, where $\boldsymbol{0}=(0,0,0)$. In our case, the normalization (2.8) gives us $f_{{\boldsymbol{0}}}=1$. In a similar manner, we can see from the other two equations in (2.7) and (2.8) that

[TABLE]

where ${e}_{i}$ is a three dimensional index whose $i$-th entry equals $1$ and other entries equal zero. The heat flux and stress tensor are related to the coefficients by

[TABLE]

3 Approximation of quadratic collision term

In order to investigate the evolution of the coefficients $f_{{\alpha}}$ in the expansion (2.11), we shall expand the collision term under the same function space. The expansion of collision operator of the linear type is rather straightforward. As an example, the explicit form of expansion of (2.10) in three dimensional case is

[TABLE]

which comes as a consequence of the property that Hermite polynomials can diagonalize the linear FP operator. It is also intrinsically implied by the fact that Fokker-Planck equation can be used in the context of stochastic process while Hermite polynomials play a crucial role in Brownian motion, but we shall not take the stochastic perspective here.

We shall first discuss the series expansion of the quadratic collision term $\mathcal{Q}[f]$ defined in (2.2), and then combine the quadratic result with the linear-type collision operators to construct collision models with better accuracy and less computational complexity.

3.1 Series expansions of quadratic collision terms

Suppose the binary collision term $\mathcal{Q}[f]$ is expanded into the following form

[TABLE]

Due to the orthogonality of Hermite polynomials, we get

[TABLE]

where the last equality can be derived by inserting (2.11) into (2.2), and

[TABLE]

The above formula is of an extremely complex form, with the evaluation of every single coefficient requiring a six dimensional integration, as well as differential operations. Granted this can be computed by numerical quadrature; the computational cost would be unbearble for getting all these coefficients. Recently, in [27], a strategy to simplify the above integral is introduced for the Coulombian case $\gamma=-3$, and the explicit values are given with small indices. In order to deal with this integral, we give the explicit expressions of all the coefficients $A_{{\alpha}}^{{\lambda},{\kappa}}$ and enlarge the applicable region of these expressions to $\gamma>-5$ for the quadratic collision kernel, which incorporates the domain of definition for $\gamma$ in the IPL model. The main results are summarized in the following theorem:

Theorem 1**.**

The expansion coefficients of the collision operator $\mathcal{Q}[f](\boldsymbol{v})$ defined in (3.3) have the following form:

[TABLE]

where $p=(p_{1},p_{2},p_{3})^{T}$ is a three-dimensional multi-index and

[TABLE]

the sum is taken for the indices in the range if and only if each subindex is non-negative.

The coefficients $a_{pq}^{\lambda\kappa}$ and $B_{p}^{q}(\gamma,s,t)$ are defined by

[TABLE]

and

[TABLE]

where

[TABLE]

The proof of Theorem 1 can be found in Appendix 7.1. Hence, we only have to compute (3.9). When $\gamma>-3$, it can be computed directly by the recursive formula of the Hermite Polynomials following the method in [35]. However, for the Coulombian case $\gamma=-3$, the recursive formula can not be adopted directly due to the singularity induced by the small value of $\gamma$. In [27], the Coulombian case $\gamma=-3$ is evaluated by adopting the special form of the quadratic collision term there. In the next section, we will introduce a new method to deal with the super singularity for a large region of $\gamma$.

3.2 Derivation of exact coefficients in super singular integral

In this section, we will introduce a different method to calculate these coefficients exactly and the applicable area of $\gamma$ is enlarged as well. In order to deal with the singularity , Burnett polynomials, products of Sonine polynomials and solid spherical harmonics [20], are utilized here. Burnett polynomials are introduced in [9] to approximate the distribution function of Boltzmann equation, and was adopted in [10, 18] to reduce the quadratic collision operator. To be concrete, the normalized form of the Burnett polynomials is

[TABLE]

where the index $\hat{{\alpha}}$ is defined as

[TABLE]

Here $L_{n}^{(\beta)}(x)$ is the Laguerre polynomials

[TABLE]

and $Y_{l}^{m}(\boldsymbol{n})$ is spherical harmonics

[TABLE]

with $P_{l}^{m}$ the associate Legendre polynomial

[TABLE]

By the orthogonality of Laguerre polynomials and spherical harmonics, one can find that

[TABLE]

In order to reduce complexity, the symmetry of $G_{st}(\gamma,p,q)$, which is stated in Lemma 2, is utilized first to reduce the cost of computation and storage.

Lemma 2**.**

For the expressions $G_{st}(\gamma,p,q)$, it holds that

[TABLE]

and

[TABLE]

Here $\Pi_{i}^{j}p$ is a permutation operator which exchanges the $i$-th and $j$-th entries of $p$.

Based on Lemma 2, we only have to compute two cases $G_{33}(\gamma,p,q)$ and $G_{13}(\gamma,p,q)$. In order to handle the singularity in $G_{st}(\gamma,p,q)$, Hermite polynomials in (3.9) is expressed by a linear combination of the Burnett polynomials, precisely

[TABLE]

where $|\hat{{\alpha}}|_{B}=\hat{\alpha}_{1}+2\hat{\alpha}_{3}.$ Since both Hermite and Burnett polynomials are orthogonal polynomials associated with the same weight function, thus the coefficients $C_{\hat{{\alpha}}}^{{\alpha}}$ defined in (3.13) are nonzero only when the degrees of $H^{{\alpha}}$ and $B_{\hat{{\alpha}}}$ are equal, precisely $|\hat{{\alpha}}|_{B}=|{\alpha}|$. The detailed algorithm to compute the coefficient $C_{\hat{{\alpha}}}^{{\alpha}}$ can be found in [10] and we also explain that briefly in Appendix 7.3. With the help of the Burnett polynomials, we can finally get the exact value of $G_{st}(\gamma,p,q).$

Proposition 3**.**

When $\gamma>-5$, $G_{st}(\gamma,p,q)$ defined in (3.9) can be simplified as

[TABLE]

where

[TABLE]

where $n_{s}$ and $n_{t}$ are the $s$-th and $t$-th entries of the unit vector $\boldsymbol{n}$ in spherical coordinates $\boldsymbol{n}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)^{T}$. The parameters in (3.14) are defined as

[TABLE]

and

[TABLE]

Proof.

Substituting (3.13) into (3.9), and adopting the formula introduced in [30, eq.(10)]

[TABLE]

we have thus validated this proposition. ∎

Finally, as stated previously in Lemma 2, we only need to compute $G_{33}(\gamma,p,q)$ and $G_{13}(\gamma,p,q)$. Therefore, only these two corresponding cases of (3.15) are discussed in the following theorem and the proof is presented in Appendix 7.2.

Theorem 4**.**

Define $\eta_{lm}^{\mu}$ as

[TABLE]

Then the coefficients $F_{13}(l_{1},m_{1},l_{2},m_{2})$ and $F_{33}(l_{1},m_{1},l_{2},m_{2})$ have the following explicit form

[TABLE]

The above analysis shows that for the FPL collision operator, the coefficients $A_{{\alpha}}^{{\lambda},{\kappa}}$ can be calculated exactly for $\gamma>-5$, which makes it much easier to build the high order scheme to numerically solve FPL equation. Moreover, this algorithm for the coefficients here is readily applicable for offline numerical evaluation, the effectiveness of which is corroborated by our numerical examples.

4 Construction of novel collision model

Until now, we already obtain a complete algorithm to calculate the coefficients $A_{{\alpha}}^{{\lambda},{\kappa}}$, which can be utilized either to discretize the quadratic collision term or to construct new collision models. We will now discuss both topics.

4.1 Discretization of homogeneous FPL equation

We will discrete the homogeneous FPL equation by the Galerkin spectral method in terms of the expansion of the distribution function (2.11). For any positive integer $M$, we define as the functional space of numerical solution

[TABLE]

where $I_{M}=\{(\alpha_{1},\alpha_{2},\alpha_{3})|0\leqslant|{\alpha}|\leqslant M,\alpha_{i}\in\mathbb{N},i=1,2,3\}$. Then the semi-discrete function $f_{M}(t,\cdot)\in\mathcal{F}_{M}$ satisfies

[TABLE]

Suppose

[TABLE]

The equations (3.2) and (3.3) imply that the variational form (4.2) is equivalent to the ODE system below:

[TABLE]

We thus obtain the formulation of the ODE system in its full form (4.4), with the help of the exact coefficients $A_{{\alpha}}^{{\lambda},{\kappa}}$ for all ${\alpha},{\lambda},{\kappa}\in I_{M}$. With $M$ fixed, these coefficients need to be computed only once, and then can be used repeatedly for multiple numerical examples.

4.2 Approximation to general collision model

In practical computation, the storage cost of computing the coefficients $A_{{\alpha}}^{{\lambda},{\kappa}}$ is formidably expensive, as the number of coefficients increases significantly with $M$ increasing. Moreover, the computational cost $O(M^{9})$ is an issue especially when solving the spatially non-homogeneous problems.

To overcome this difficulty, the method in [35] is utilized to reduce the computational cost, precisely that the coefficients $A_{{\alpha}}^{{\lambda},{\kappa}}$ for a small number $M_{0}$ are computed and stored, when the computational cost for solving (4.4) is acceptable. As for ${\alpha}\not\in I_{M_{0}}$, we apply the linear model (2.10) brought up by Villani and compute as:

[TABLE]

Combining (4.4) and (4.5), we obtain a novel collision operator

[TABLE]

where $P_{M_{0}}$ is the orthogonal projection from $\mathcal{F}$ onto $\mathcal{F}_{M_{0}}$. After applying spectral method to this collision operator in the functional space $\mathcal{F}_{M}$, where $M$ is chosen to be larger than ${M_{0}}$ the final ODE system for the new model is

[TABLE]

where

[TABLE]

By now, we have obtained a series of new collision models (4.6). It can be expected that such combination could reduce the time cost significantly due to the simple form of the linear FP collision operator in the Hermite basis, while at the same time manage to maintain a high level of accuracy since the evolution function already captures the most crucial information in coefficients of lower order and performs a satisfactory approximation in the other coefficients. This will be observed in the numerical examples.

5 Numerical examples

In this section, we shall present several results in our numerical computation. In all of the numerical experiments, we shall adopt the newly proposed collision operator (4.6), and solve the equation

[TABLE]

numerically for some positive integer $M_{0}$. Such an equation is solved by the Galerkin spectral method for solution defined in the functional space $\mathcal{F}_{M}$, with $M$ chosen to be greater than $M_{0}$. Namely we solve the system of ODE (4.8). As for the discretization in time, we use the 4th-order Runge-Kutta method in the examples, and the time step is chosen as $\Delta t=0.01$. In the examples, we shall set $\Lambda=1.$

Finally, we would like to mention that the derivation of the expansion coefficients in the Hermite basis are exact in each case by mathematical derivation instead of numerical integration in order to achieve high accuracy.

5.1 BKW solution

For the Maxwell gas $\gamma=0$, the original FPL equation admits an exact solution with the following expression:

[TABLE]

where $\tau(t)=1-0.4\exp\left(-4t\right)$. As a good approximation of the initial distribution function, we use $M=15$ ($816$ degrees of freedom) in our simulation. For the visualization purpose, we define the marginal distribution functions (MDFs)

[TABLE]

The initial MDFs are plotted in Figure 1, in which the lines for exact functions and their numerical approximation are hardly distinguishable.

Numerical results for $t=0.01$, $0.02$ and $0.06$ are given in Figures 2 and 3, respectively for the marginal distribution function $g(t,v_{1})$ and $h(t,v_{1},v_{2})$. Here we set $M_{0}$ as $M_{0}=5$ and $15$. For $M_{0}=5$, the numerical solution provides a reasonable approximation, but still with noticeable deviations, while for $M_{0}=15$, the two solutions match perfectly in all cases.

Now we consider the time evolution of the expansion coefficients. By expanding the exact solution into Hermite series, we get the exact solution for the coefficients:

[TABLE]

From (5.1), we can find that the coefficients $f_{\alpha}$ are zero for any $t$ if $1\leqslant|\alpha|\leqslant 3$. Hence we will focus on the coefficients $f_{400}$ and $f_{220}$ here. Figure 4 gives the comparison between the numerical solution and the exact solution for these two coefficients. In both plots, all the three lines coincide perfectly.

5.2 Bi-Gaussian initial data

In this example, we perform the numerical test for the Bi-Gaussian problem. Here the Coulombian case $\gamma=-3$ are tested. The initial distribution function is

[TABLE]

In this numerical test, we use $M=20$ which gives a good approximation of the initial distribution function (see Figure 5).

For this example, we consider the three cases $M_{0}=5,10,15$, and the corresponding one-dimensional marginal distribution functions at $t=0.4$, $1$ and $3$ are given in Figure 6. In all the results, the numerical results are converging to the reference solution as $M_{0}=15$, and the lines for $M_{0}=10$ and $M_{0}=15$ are very close to each other. To get a clearer picture, similar comparisons of two-dimensional results are also provided in Figure 7.

Now we consider the evolution of stress tensor and heat flux. In this example, we always have $\sigma_{11}=-2\sigma_{22}=-2\sigma_{33}$ and $q_{i}=0,i=1,2,3$. Therefore we focus only on the evolution of $\sigma_{11}$, which is plotted in Figure 8. It can be seen that even for $M_{0}=5$, the evolution of the stress tensor is almost exact, where the distribution function is not approximated very well. The three lines are all on top of each other.

In order to test the computational capacity of our new model, the same example with a very small $\gamma$ as $\gamma=-4.9$ is tested. Here we also set $M=20$, and choose the numerical results with $M_{0}=15$ as the reference solution. Figure 9 shows the marginal distribution $h(t,v_{1},v_{2})$ at $t=0.4,1$ and $3$ with $M_{0}=5$ and $10$. It illustrates that when $\gamma$ equals $-4.9$, the numerical solutions are converging to the reference solution as $M_{0}=15$, and that the solution with $M_{0}=10$ is almost the same as the reference solution. The time evolution of $\sigma_{11}$ is plotted in Figure 10, where the three results are also on top of each other, even with $M_{0}=5$. Moreover, from Figure 9 and 10, we can find that the time evolution of the distribution function with $\gamma=-4.9$ is slower than that with $\gamma=-3$, which is also consistent with the form of the FPL collision operator.

5.3 Rosenbluth problem

In this example, the Rosenbluth problem is tested. Also the Coulombian case $\gamma=-3$ and the case $\gamma=-4.9$ are tested. The initial condition is from [25] as

[TABLE]

The parameter $A,B$ are standardized to satisfy the condition that the initial density and temperature all equal $1$, precisely

[TABLE]

where $a=\pi(3\sqrt{\pi}({\rm erf}(1)+1+2/e)$ and $b=\pi(9.5\sqrt{\pi}({\rm erf}(1)+1)+7/e)$. Here $e$ is the Euler number and $\rm erf(x)=\frac{1}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}\,\mathrm{d}t$ is the error function.

In order to approximate the initial distribution function well, in this numerical test, we set $M=40$. The initial MDFs are plotted in Figure 11, which illustrates the perfect numerical approximation to the exact distribution function.

For this example, also the three cases $M_{0}=5,10,15$ are tested. The numerical solution with $M_{0}=15$ is treated as the reference solution. The corresponding one and two dimensional marginal distribution functions for the Coulombian $\gamma=-3$ are shown in Figure 12 and 13, where the numerical solutions are converging to the reference solution and that with $M_{0}=10$ is almost the same as the reference solution.

Moreover, our new model can also approximate the $\gamma=-4.9$ case well. The two dimensional marginal distribution functions for $\gamma=-4.9$ are presented in Figure 14. The numerical solution with $M_{0}=15$ is also chosen as the reference solution. Similar to the example in Sec 5.2, the time evolution of the distribution function with $\gamma=-4.9$ is also slower than that with $\gamma=-3$. Further more, though there are distinct differences between the numerical solution $M_{0}=5$ and $M_{0}=15$, the numerical solutions are converging to the reference solutions $M_{0}=15$.

6 Conclusion

In this paper, we focus on applying the Hermite spectral method to develop an efficient and accurate way of approximating and numerically solving the Fokker-Planck-Landau equation. Basic properties of Hermite polynomials are utilized to obtain a simplified expression of the coefficients, which renders the numerical method feasible. Burnett polynomials are introduced to deal with the super singular integral in the computation. This method could cover more practical cases up to $\gamma>-5$.

A novel collision model is built with a combination of quadratic collision model and the linearized collision model brought up by C. Villani [34]. The numerical experiments validate the efficiency of this new model. With the model introduced, the numerical solutions are of a high accuracy, as well as an affordable computational cost. This method should be further validated in the numerical tests for the full FPL equation with spatial variables, which will be one of the future works.

Acknowledgements

Ruo Li is supported by the National National Scientific Foundation of China (Grant No. 91630310) and Science Challenge Project (No. TZ2016002). Yanli Wang is supported by the National Natural Scientific Foundation of China (Grant No. 11501042), and Chinese Postdoctoral Science Foundation of China (2018M631233).

7 Appendix

7.1 Proof of Theorem 1

We shall present of proof of Theorem 1 here. In order to prove Theorem 1, we first introduce the lemma below:

Lemma 5**.**

Let $\boldsymbol{v}=\boldsymbol{h}+\boldsymbol{g}/2$ and $\boldsymbol{w}=\boldsymbol{h}-\boldsymbol{g}/2$. It holds that

[TABLE]

where the coefficients $a_{\alpha^{\prime}\kappa^{\prime}}^{\alpha\kappa}$ are defined in (3.7).

Corollary 1**.**

Let $\boldsymbol{v}=\boldsymbol{h}+\boldsymbol{g}/2$. We have

[TABLE]

The proof of Lemma 5 can be found in [35]. Next we will prove Theorem 1.

Proof of Theorem 1.

Using an integration by parts and the recursion formula of Hermite polynomials

[TABLE]

the coefficients $A_{\alpha}^{\lambda,\kappa}$ (3.4) can be simplified as

[TABLE]

where $G_{st}(\boldsymbol{v})=-v_{s}v_{t}+\delta_{st}|\boldsymbol{v}|^{2}$. Further simplification of (7.2) follows the method in [35], where the velocity of the mass center is defined as $\boldsymbol{h}=(\boldsymbol{v}+\boldsymbol{v}^{\ast})/2$ and the relative velocity is defined as $\boldsymbol{g}=\boldsymbol{v}-\boldsymbol{v}^{\ast}$. Hence, it holds that

[TABLE]

Combing Lemma 5, (7.2) and (7.3), the integral in $A_{\alpha}^{\lambda,\kappa}$ can be rewritten with respect to $\boldsymbol{g}$ and $\boldsymbol{h}$

[TABLE]

Using the orthogonality of Hermite polynomials, we can finally prove Theorem 1.

∎

7.2 Proof of Theorem 4

In order to prove Theorem 4, we will introduce the lemma below:

Lemma 6**.**

For three spherical harmonics $Y_{l}^{m}$, $Y_{l_{1}}^{m_{1}}$ and $Y_{l_{2}}^{m_{2}}$, if $m\neq m_{1}+m_{2}$, or $l\not\in[|l_{1}-l_{2}|,l_{1}+l_{2}]$, then

[TABLE]

Especially, we have

[TABLE]

where $\eta_{lm}^{\mu}$ is defined in (3.18).

The result of this lemma can be found in Section 12.9 of [1].

Proof of Theorem 4.

Noting that

[TABLE]

Based on Lemma 6 and the property of spherical harmonic $\overline{Y_{l}^{m}(\boldsymbol{n})}=(-1)^{m}Y_{l}^{-m}(\boldsymbol{n})$, we can derive the results in Theorem 4 with the orthogonality property of spherical harmonics

[TABLE]

∎

7.3 Computation of Coefficients

$C_{\hat{{\alpha}}}^{{\alpha}}$

In this section, we will briefly introduce the algorithm to calculate $C_{\hat{{\alpha}}}^{{\alpha}}$, and the original algorithm is in [10].

Define

[TABLE]

and the recursive formula of the basis functions [11] is

[TABLE]

where $\eta_{lm}^{\mu}$ is defined in (3.9) and we set $B_{\hat{{\alpha}}}(\boldsymbol{v})=0$ if $|\hat{\alpha}_{2}|>\hat{\alpha}_{1}$ or either of $\hat{\alpha}_{1},\hat{\alpha}_{3}$ is negative. Based on the recursion formula of Hermite polynomials

[TABLE]

we can get the recursive formula to compute $C_{\hat{{\alpha}}}^{{\alpha}}$, precisely

[TABLE]

where $|{\alpha}|=|\hat{{\alpha}}|_{B}$ and

[TABLE]

As is stated in [10], we solve all the coefficients $C_{\hat{{\alpha}}}^{{\alpha}}$ by the order of $|{\alpha}|$, so that the right-hand sides of (7.12) are always known. The initial condition and the boundary conditions are $C_{\boldsymbol{0}}^{\boldsymbol{0}}=1$ and $C_{\hat{{\alpha}}}^{{\alpha}}=0$ if $|\hat{\alpha}_{2}|>\hat{\alpha}_{1}$ or either of $\hat{\alpha}_{1},\hat{\alpha}_{3}$ is negative. Moreover, the time complexity for computing all the coefficients $C_{\hat{{\alpha}}}^{{\alpha}}$ with $|\hat{{\alpha}}|_{B}=|{\alpha}|\leqslant M$ is $O(M^{5})$.