Solving Vlasov Equations Using NR Method
Abstract
In this paper, we propose a moment method to numerically solve the Vlasov equations using the framework of the NR method developed in [6, 8, 7] for the Boltzmann equation. Due to the same convection term of the Boltzmann equation and the Vlasov equation, it is very convenient to use the moment expansion in the NR method to approximate the distribution function in the Vlasov equations. The moment closure recently presented in [5] is applied to achieve the globally hyperbolicity so that the local wellposedness of the moment system is attained. This makes our simulations using high order moment expansion accessible in the case of the distribution far away from the equilibrium which appears very often in the solution of the Vlasov equations. With the moment expansion of the distribution function, the acceleration in the velocity space results in an ordinary differential system of the macroscopic velocity, thus is easy to be handled. The numerical method we developed can keep both the mass and the momentum conserved. We carry out the simulations of both the VlasovPoisson equations and the VlasovPoissonBGK equations to study the linear Landau damping. The numerical convergence is exhibited in terms of the moment number and the spatial grid size, respectively. The variation of discretized energy as well as the dependence of the recurrence time on moment order is investigated. The linear Landau damping is well captured for different wave numbers and collision frequencies. We find that the Landau damping rate linearly and monotonically converges in the spatial grid size. The results are in perfect agreement with the theoretic data in the collisionless case.
Keywords: Vlasov equations; NR method; Landau damping
1 Introduction
The Vlasov equation is a differential equation describing the time evolution of the distribution function of plasma consisting of charged particles with the longrange (for example, Coulomb) interaction. The equation was first suggested for the description of plasma by A. Vlasov in 1938 [16]. Due to the presence of longrange Coulomb interaction in the plasma, Vlasov started from the Vlasov equation, which is a collisionless Boltzmann equation, and used a selfconsistent collective field created by the charged plasma particles to get the VlasovMaxwell equations [11]. The VlasovPoisson (VP) equations are an approximation of the VlasovMaxwell equations in the nonrelativistic zeromagnetic field limit. The VP equations are used to describe various phenomena in plasma, in particular Landau damping and the distributions in a double layer plasma. The VlasovPissionBGK (VB) equations are also studied with a collision term presented as the BGK term. They are the simplest kinetic equations which correctly describe the essential features of collective and dissipative (entropyproducing) particle interactions in semiconductor plasmas [1, 2, 17, 18].
Due to the complex phenomena in the plasma, numerical simulation plays an important role in the study of the Vlasov and related equations. There are several kinds of methods to solve the Vlasov equations. The finite element method was proposed in [19, 20], which can be used to handle complicated boundary problems but inconvenient to solve the high dimension equations [10]. Meanwhile, the particleincell (PIC) method [3] which used a finite number of macroparticles to approximate the plasma is easy to be implemented, and the method to discretize the Vlasov equations on a mesh of phase space was introduced to remedy the problem in PIC that the inherent artificial discrete particle noise made the description inaccuracy. The semiLagrangian method [15] and the cubic interpolated propagation method [12] were also used to solve the Vlasov equations. However, the first method destroyed the local characters due to the reconstruction and the second one was quite expensive for the storage of the distribution function and its derivatives. In [10], Filbet put forward a new method to deal with the force term, which made the scheme keep the mass and energy conserved. Recently, an approach based on the moment method has been proposed in [13], and therein the distribution function was expanded using the Hermite polynomials with a prescribed macroscopic velocity chosen as the expansion center and a prescribed thermal velocity as the scaling factor at different locations.
In the past years, a regularized moment method was developed in [6] to numerically solve the Boltzmann equation. This method adopts the Hermite polynomial expansion to approximate the distribution function, with the basis function shifted by the local macroscopic velocity and scaled by the square root of the local temperature. The approximated distribution function is used to directly solve the Boltzmann equation without the deduction of the moment system up to arbitrary order of moments. The method therein was further explored as the NR method [8, 7] by introducing the regularization term using asymptotic expansion in term of the mean free path. Recently, a new regularization method [5, 4] was derived with the guarantee that the regularized moment system is globally hyperbolic. Due to the locally wellposedness provided by the global hyperbolicity, it is eventually accessible that approximating the distribution function far away from the equilibrium distribution by the stable simulation using large number of moments. Inspired by this progress, we in this paper develop the NR method to study the Vlasov equations, which is similar to the Boltzmann equation in the convection term, while the distribution function is much farther away from the equilibrium state than the gas flows, due to the longrange Coulomb interactions.
Here we are focusing on the VP and VB equations. With the moment expansion in the Vlasov equations, the convection part is smoothly handled by the original NR method for the Boltzmann equation. We discretize the regularized term given in [4] directly using the central difference scheme. Our deduction shows that the electric potential brings us an acceleration on the macroscopic velocity, thus it is very convenient to be numerically integrated. Currently, the collision term under our consideration is the simple BGK model to avoid distraction. The approximated electric potential is obtained by a threepoint central scheme, and the numerical method keeps both the mass and the momentum conserved. By the numerical resolution study, it is exhibited that our method is numerically converged by the comparison of Landau damping rates obtained using different spatial grid size and number of moments. With the increasing of the number of moments, the recurrence time is almost linearly related to the square root of the moment expansion order. The discretized energy of our method only varies slightly in comparison with the overall energy. Since the Landau damping rate is affected by the wave number and the collision frequency, different wave numbers and collision frequencies are extensively studied. The wave numbers ranging from 0.2 to 0.5 are simulated, and the collision frequencies are taken as 0.0 (the collisionless case), 0.01 and 0.05. The results show that our method can capture the linear Landau damping very well and the damping rates obtained is in quantitative agreement with the theoretic data. We find with surprise that the numerical Landau damping rate is linearly and monotonically converged in term of the spatial grid size. Particularly, the numerical damping rate converges to the theoretic data perfectly in the collisionless case. This observation inspires us to predict the Landau damping rates by the extrapolation of our numerical damping rates with presence of the collision term.
The layout of this paper is as follows: in Section 2, the regularized moment system is deduced for the VlasovBGK equations. In Section 3, we present the detailed procedure of the numerical method. In Section 4, the numerical examples including the numerical resolution study and the linear Landau damping simulation with different parameters are presented. Some concluding remarks are given in the last section.
2 Regularized Moment System
Let , which depends on time , position and the microscopic velocity , be the distribution function depicting the motion of particles. It is governed by the VB equations
(2.1) 
where denotes the collision frequency, and is the electric force produced by the selfconsistent electric filed :
(2.2) 
where is the electric potential produced by the particles; , , and stand for the density, the single charge, the mass of one particle and the electric constant respectively; is the Maxwellian defined as
(2.3) 
where the parameter is the thermal velocity [13], is the macroscopic velocity and is the same as that in (2.2). is related to by
(2.4) 
In the case of , we get the VP equations. The relations between the macroscopic variables and the distribution function are deduced as
(2.5)  
(2.6)  
(2.7) 
The conservation of mass, momentum and total energy are all valid for the VB and VP equations. Multiplying the equation (2.1) by and , direct integration with and gives us
(2.8)  
(2.9) 
Multiplying the equation (2.1) by and integrating by parts, we get the conservation of the total energy for the system (2.1) and (2.2):
(2.10) 
2.1 Hermite expansion of the distribution function
Following the method in [6, 9], we expand the distribution function into Hermite series as
(2.11) 
where is a threedimensional multiindex, and
(2.12) 
The basis functions are defined as
(2.13) 
where is the Hermite polynomial
(2.14) 
For convenience, is taken as zero if , thus is zero when any component of is negative. With such an expansion, the Maxwellian can be written as
(2.15) 
and the BGK collision term is written as
(2.16) 
The definition of shows that each basis function is an exponentially decreasing function multiplied by a multidimensional Hermite polynomial shifted by the local macroscopic velocity and scaled by the square root of the local thermal velocity . For any vector and positive number , if the distribution function is expanded as
(2.17) 
then the following relations hold
(2.18a)  
(2.18b)  
(2.18c) 
In the case of and in (2.18), the following relations between the coefficients can be verified:
(2.19) 
Moreover, direct calculations give us the relations between the coefficients in (2.11) as
(2.20)  
(2.21) 
where , and are related to by
(2.22)  
(2.23) 
2.2 Moment expansion of Vlasov equation
To get the moment system, we substitute the expansion (2.11) into (2.1) and then match the coefficients for the same basis functions. Taking the temporal and spatial derivatives directly on the basis functions , the term with the expansion (2.11)
is expanded as
(2.24) 
the acceleration term is expanded as
(2.25) 
and the collision term is expanded as
(2.26) 
Substituting these expansions back into (2.1), and collecting coefficients for the same basis functions, we get the following general moment equations with a slight rearrangement:
(2.27) 
where is defined by
(2.28) 
Following the method in [7], we deduce the mass conservation in the case of as
(2.29) 
If we set , with , (2.27) reduces to
(2.30) 
which is simplified as
(2.31) 
Then we consider the case of . Multiplying on both sides of (2.1), and integrating with respect to on , we have
(2.32) 
Remark 1.
2.3 Closure of the moment system
With a truncation of (2.11), (2.34) will result in a finite moment system. Precisely, we let be a positive integer and only the coefficients in the set are considered. Let denotes the linear space spanned by all ’s with , and the expansion (2.11) is truncated as
(2.35) 
with and . The moment equations which contain with are disregarded in (2.34). Then, (2.29), (2.31) and (2.34) with lead to a system with finite number of equations.
Due to the presence of the terms with , , the moment system we obtained is not closed yet. We rewrite (2.34) into the form below:
(2.36) 
where in the case of ,
and in the case of ,
Clearly, the moments in term with are not in the set of moments , and have to be substituted by some expressions consisting of lower order moments to make the moment system closed. If we simply let , the Gradtype system associated with the moment set is obtained. It is known that the Gradtype system is not locally wellposed due to the lack of the global hyperbolicity, which results in numerical blowup when the distribution function is far away from the equilibrium state. The regularization method in [4, 5] is adopted here to achieve a globally hyperbolic moment closure. The st order terms are substituted as bellow
(2.37) 
Let denote the regularization term based on the characteristic speed correction in [4, 5] for
(2.38) 
Substituting the st order term with the regularization term , the moment equations are revised as
(2.39) 
with for . If the distribution function only depends on in the spatial direction, we have that
(2.40) 
Remark 2.
The regularized moment system (2.39) is not able to be written into a conservation law, for the presence of the regularization term. If we let with , the system changes into the conservative Gradtype system.
3 Numerical Method
The numerical scheme for the regularized moment system (2.39) is a natural extension of the method in [6]. By a standard fraction step method, we split the VB equations into the following parts:

the convection step:
(3.1) 
the acceleration step:
(3.2) (3.3) 
the collision step:
(3.4)
We observe that only (2.31) contains the electric force in the governing equations. Thus the governing equations of the acceleration part turn into
(3.5)  
(3.6) 
Here we restrict our study in 1D spatial space. The numerical scheme adopted in the direction is the standard finite volume discretization. Suppose to be a uniform mesh in , and each cell is identified by an index . For a fixed and ,
(3.7) 
The numerical solution which is the approximation of the distribution function at is denoted as
(3.8) 
where is the approximation over the cell at the th time step and has the following Hermite expansion form as
3.1 The convection step
The equation (3.1) is discretized as
(3.9) 
where is the contribution of the term in (2.39) without considering the acceleration, and is the contribution of the term in (2.39). Noticing that the term results in the conservative part in the Gradtype moment system, its contribution is discretized in the conservative formation as
(3.10) 
where is the numerical flux between cell and at . We use the HLL scheme in our numerical experiments following [8], which reads:
(3.11) 
where and are the fastest signal speeds [5] as
(3.12) 
where is the greatest zero of , is the first component of the macroscopic velocity , and is the thermal velocity. The formula of the signal speed is also used to determine the time step by the CFL condition. Two points remain unclear in the numerical flux. The first one is how to calculate . This is managed according to the recursion relation of Hermite polynomials:
(3.13) 
where , is the first entry of , and are the mean macroscopic velocity and thermal velocity in the th cell. Since when , no longer exists in the space . By simply dropping the terms with , we project back into , since when , is orthogonal to with respect to the inner product
(3.14) 
The other point is how to add up two distribution functions lying in and respectively. The proposition in [6] is referred to solve it.
Proposition 1.
Suppose can be represented by
(3.15) 
For some and , satisfies
(3.16) 
where and are given below. And , .
(3.17) 
Let
(3.18) 
Then and satisfies
(3.19) 
The proposition provides an algorithm to project a function in to the space . Assuming and , we can find an as a representation of in the sense of (3.19). Thus, to add up and , we first find an as an approximation of in the sense of (3.19), and then add up the coefficients of and for the same basis function.
The regularization term appears only in the moment equations containing with . Therefore, only when , we have to calculate . We simply use the central difference scheme to approximate the spatial derivatives in (2.38):
(3.20) 
Then we get the numerical approximation for with , as
(3.21) 
3.2 The acceleration and collision step
The acceleration step is performed by solving
(3.22) 
For the time step , (3.22) is approximated as
(3.23) 
where is the first entry of the macroscopic velocity in the th cell after the convection step at . And is the first entry of the electric force in the th cell after the convection step at . In the 1D spatial space, the Poisson equation (2.2) reduces into a second order ODE as
(3.24) 
The threepoint central difference scheme is used to discretize the potential equation
(3.25) 
where is the density in the th cell after the convection step, noticing that the density is not updated in the acceleration step and the collision step. The central difference is used to approximate and