Solving Vlasov Equations Using NRxx Method

# Solving Vlasov Equations Using NRxx Method

Zhenning Cai,   Ruo Li   and Yanli Wang School of Mathematical Sciences, Peking University, Beijing, China, email: caizn@pku.edu.cn.HEDPS & CAPT, LMAM & School of Mathematical Sciences, Peking University, Beijing, China, email: rli@math.pku.edu.cn.CAPT, Beijing International Center for Mathematical Research, & School of Mathematical Sciences, Peking University, Beijing, China, email: wangyanliwyl@gmail.com.
###### 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 well-posedness 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 Vlasov-Poisson equations and the Vlasov-Poisson-BGK 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 long-range (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 long-range Coulomb interaction in the plasma, Vlasov started from the Vlasov equation, which is a collisionless Boltzmann equation, and used a self-consistent collective field created by the charged plasma particles to get the Vlasov-Maxwell equations [11]. The Vlasov-Poisson (V-P) equations are an approximation of the Vlasov-Maxwell equations in the nonrelativistic zero-magnetic field limit. The V-P equations are used to describe various phenomena in plasma, in particular Landau damping and the distributions in a double layer plasma. The Vlasov-Pission-BGK (V-B) 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 (entropy-producing) 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 particle-in-cell (PIC) method [3] which used a finite number of macro-particles 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 semi-Lagrangian 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 well-posedness 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 long-range Coulomb interactions.

Here we are focusing on the V-P and V-B 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 three-point 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 Vlasov-BGK 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 V-B equations

 ∂f∂t+v⋅∇xf+F(t,x,v)⋅∇vf=ν(fM−f), (2.1)

where denotes the collision frequency, and is the electric force produced by the self-consistent electric filed :

 F(t,x,v)=qmE(t,x),E(t,x)=−∇x ϕ(t,x),−Δxϕ=qρϵ0, (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

 fM=ρ(t,x)(2πuth(t,x))3/2exp(−|v−u(t,x)|22uth(t,x)), (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

 ∫R3f(v)⎛⎜⎝1v|v|2⎞⎟⎠dv=∫R3fM(v)⎛⎜⎝1v|v|2⎞⎟⎠dv. (2.4)

In the case of , we get the V-P equations. The relations between the macroscopic variables and the distribution function are deduced as

 ρ(t,x)=∫R3f(t,x,v)dv, (2.5) ρ(t,x)u(t,x)=∫R3vf(t,x,v)dv, (2.6) ρ(t,x)|u(t,x)|2+3ρ(t,x)uth(t,x)=∫R3|v|2f(t,x,v)dv. (2.7)

The conservation of mass, momentum and total energy are all valid for the V-B and V-P equations. Multiplying the equation (2.1) by and , direct integration with and gives us

 ddt∫R3×R3f(t,x,v)dxdv=0,t∈R+, (2.8) ddt∫R3×R3vf(t,x,v)dxdv=0,t∈R+. (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):

 ddt(∫R3×R3f(t,x,v)|v|2dxdv+∫R3|E(t,x)|2dx)=0,t∈R+. (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

 f(v)=∑α∈N3fαHuth,α(ξ), (2.11)

where is a three-dimensional multi-index, and

 ξ=v−u√uth. (2.12)

The basis functions are defined as

 Huth,α(ξ)=3∏d=11√2πu−αd+12thHeαd(ξd)exp(−ξ2d2), (2.13)

where is the Hermite polynomial

 Hen(x)=(−1)nexp(x22)dndxnexp(−x22). (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

 fM(v)=ρHuth,0(ξ), (2.15)

and the BGK collision term is written as

 ν(fM−f)=−ν∑|α|≥1fαHuth,α(ξ). (2.16)

The definition of shows that each basis function is an exponentially decreasing function multiplied by a multi-dimensional 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

 f(v)=∑α∈N3f′αHu′th,α(ξ′),ξ′=v−u′√u′th, (2.17)

then the following relations hold

 ρ =f′0, (2.18a) ρu =ρu′+(f′ed)Td=1,2,3, (2.18b) ρ|u|2+3ρuth =2ρu⋅u′−ρ|u′|2+3∑d=1(u′thf′0+2f′2ed). (2.18c)

In the case of and in (2.18), the following relations between the coefficients can be verified:

 f0=ρ(t,x),fei=0,3∑d=1f2ed=0,i=1,2,3. (2.19)

Moreover, direct calculations give us the relations between the coefficients in (2.11) as

 qi=2f3ei+3∑d=1f2ed+ei, (2.20) pij−13δij3∑d=1pdd=(1+δij)fei+ej. (2.21)

where , and are related to by

 qi=12∫R3|v−u|2(vi−ui)fdv, (2.22) pi,j=∫R3(vi−ui)(vj−uj)fdv. (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)

 ∂f∂t+v⋅∇xf

is expanded as

 ∑α∈N3{(∂fα∂t+3∑d=1∂ud∂tfα−ed+12∂uth∂t3∑d=1fα−2ed)+3∑j=1[(uth∂fα−ej∂xj+uj∂fα∂xj+(αj+1)∂fα+ej∂xj)+3∑d=1∂ud∂xj(uthfα−ed−ej+ujfα−ed+(αj+1)fα−ed+ej)+12∂uth∂xj3∑d=1(uthfα−2ed−ej+ujfα−2ed+(αj+1)fα−2ed+ej)]}Huth,α(ξ−u√uth), (2.24)

the acceleration term is expanded as

 −∑α∈N33∑d=1Fdfα−edHuth,α(ξ−u√uth), (2.25)

and the collision term is expanded as

 ν∑|α|>1fαHuth,α(ξ−u√uth). (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

 δ(α)={0,if  |α|⩾2,1,otherwise. (2.28)

Following the method in [7], we deduce the mass conservation in the case of as

 ∂f0∂xj+3∑j=1(uj∂f0∂xj+f0∂uj∂xj)=0. (2.29)

If we set , with , (2.27) reduces to

 f0(∂ud∂t+3∑j=1uj∂ud∂xj−Fd)+f0∂uth∂xd+uth∂f0∂xd+3∑j=1(δjd+1)∂fed+ej∂xj=0, (2.30)

which is simplified as

 f0(∂ud∂t+3∑j=1uj∂ud∂xj−Fd)+3∑j=1∂pjd∂xj=0. (2.31)

Then we consider the case of . Multiplying on both sides of (2.1), and integrating with respect to on , we have

 f0(∂uth∂t+3∑j=1uj∂uth∂xj)+233∑j=1(∂qj∂xj+3∑d=1pjd∂ud∂xj)=0. (2.32)
###### Remark 1.

Since

 ∫R3|v−u|2∂f∂vidv=−2∫R3(vi−ui)fdv=0,i=1,2,3, (2.33)

the acceleration term does not appear in (2.32).

Substituting (2.31) and (2.32) into (2.27), we eliminate the temporal derivatives of and . Then the quasi-linear form of the moment system reads:

 (2.34)

We collect the equations (2.29), (2.31), (2.32) and (2.34) together to obtain a moment system with infinite number of equations.

### 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

 f(v)≈∑|α|⩽MfαHuth,α(v−u√uth), (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:

 ∂fα∂t+Aα+Bα=ν(1−δ(α))fα, (2.36)

where in the case of ,

 Aα=−1f03∑d=13∑j=1∂pjd∂xjfα−ed+⋯+3∑j=1(uth∂fα−ej∂xj+uj∂fα∂xj)+3∑j=1(αj+1)∂fα+ej∂xj,Bα=0,

and in the case of ,

 Aα=−1f03∑d=13∑j=1∂pjd∂xjfα−ed+⋯+3∑j=1(uth∂fα−ej∂xj+uj∂fα∂xj),Bα=3∑j=1Bα,j,Bα,j=(αj+1)∂fα+ej∂xj.

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 Grad-type system associated with the moment set is obtained. It is known that the Grad-type system is not locally well-posed due to the lack of the global hyperbolicity, which results in numerical blow-up 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

 ∂fα+ej∂xj⟶−(3∑d=1fα−ed+ej∂ud∂xj+12(D∑d=1fα−2ed+ej∂uth∂xj)),|α|=M,  j=1,2,3. (2.37)

Let denote the regularization term based on the characteristic speed correction in [4, 5] for

 ^Bα=−3∑d=1^Bα,j,^Bα,j=−(3∑d=1fα−ed+ej∂ud∂xj+12D∑d=1fα−2ed+ej∂uth∂xj). (2.38)

Substituting the -st order term with the regularization term , the moment equations are revised as

 ∂fα∂t+Aα+^Bα=ν(1−δ(α))fα, (2.39)

with for . If the distribution function only depends on in the spatial direction, we have that

 ^Bα,j=0,for  j=2,3. (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 Grad-type 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 V-B equations into the following parts:

• the convection step:

 ∂f∂t+v⋅∇xf=0, (3.1)
• the acceleration step:

 ∂f∂t+F(t,x,v)⋅∇vf=0, (3.2) F(t,x,v)=qmE(t,x),E(t,x)=−∇x ϕ(t,x),−Δxϕ=qρϵ0. (3.3)
• the collision step:

 ∂f∂t=ν(fM−f). (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

 ∂tu=F, (3.5) F(t,x,v)=qmE(t,x),E(t,x)=−∇x ϕ(t,x),−Δxϕ=qρϵ0. (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 ,

 Γh={Tj=x0+(jΔx, (j+1)Δx):j∈Z}. (3.7)

The numerical solution which is the approximation of the distribution function at is denoted as

 fnh(x,v)=fnj(v),x∈Tj, (3.8)

where is the approximation over the cell at the -th time step and has the following Hermite expansion form as

 fnj(v)=∑|α|⩽Mfnα,jHunth,j,α⎛⎜ ⎜⎝v−unj√unth,j⎞⎟ ⎟⎠.

### 3.1 The convection step

The equation (3.1) is discretized as

 fn+1,∗j(v)=fnj(v)+Kn1,j(v)+Kn2,j(v), (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 Grad-type moment system, its contribution is discretized in the conservative formation as

 Kn1,j(v)=−ΔtnΔx[Fnj+12(v)−Fnj−12(v)], (3.10)

where is the numerical flux between cell and at . We use the HLL scheme in our numerical experiments following [8], which reads:

 Fnj+12(v)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩v1fnj(v),0⩽λLj+12,λRj+12v1fnj(v)−λLj+12v1fnj+1(v)+λLj+12λRj+12[fnj+1(v)−fnj(v)]λRj+12−λLj+12,λLj+12<0<λRj+12,v1fnj+1(v),0⩾λRj+12, (3.11)

where and are the fastest signal speeds [5] as

 λLj+12=min{un1,j−CM+1√unth,j,un1,j+1−CM+1√unth,j+1},λRj+12=max{un1,j+CM+1√unth,j,un1,j+1+CM+1√unth,j+1}, (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:

 v1fnj(v)=[(unth,j)1/2(ξn1,j)+(un1,j)]∑|α|⩽Mfnj,αHnj,α(ξnj)=∑|α|≤Mfnj,α[unth,jHnj,α+e1(ξnj)+(un1,j)Hnj,α(ξnj)+α1Hnj,α−e1(ξnj)], (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

 (f,g)=∫R3f(ξ)g(ξ)exp(|ξ|22)dξ. (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

 f(v)=∑|α|⩽Mf1,αHuth,1,α(ξ1),ξ1=(v−u1)/√uth,1. (3.15)

For some and , satisfies

 ⎧⎪ ⎪⎨⎪ ⎪⎩dFαdτ=3∑d=1S2[uth,1RFα−2ed+wd√uth,1Fα−ed],∀τ∈[0,1],Fα(0)=f1,α, (3.16)

where and are given below. And , .

 R(τ)=^uth−1(^uth−1)τ+1,S(τ)=1−τR(τ)=1(^uth−1)τ+1. (3.17)

Let

 g(v)=∑|α|⩽MFα(1)Huth,2,α(ξ2),ξ2=(v−u2)/√uth,2. (3.18)

Then and satisfies

 ∫R3p(v)f(v)dv=∫R3p(v)g(v)dv,∀p(v)∈PM(R3). (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):

 ∂ud∂x≈∇hxund,j≜und,j+1−und,j−12Δx,∂uth∂x≈∇hxunth,j≜unth,j+1−unth,j−12Δx. (3.20)

Then we get the numerical approximation for with , as

 Kn2,j=−Δt∑|α|=M(α1+1)3∑d=1(fnα−ed+e1∇hxund,j+fnα−2ed+e12∇hxunth,j)Hunth,j,α⎛⎜ ⎜⎝v−unj√unth,j⎞⎟ ⎟⎠. (3.21)

### 3.2 The acceleration and collision step

The acceleration step is performed by solving

 ∂u1∂t−F1=0. (3.22)

For the time step , (3.22) is approximated as

 un+11,j=un+1,∗1,j+ΔtFn+11,j, (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

 F1(t,x,v)=qmE1(t,x),E1(t,x)=−∂ψ(t,x)∂x,−∂xxψ=qρϵ0. (3.24)

The three-point central difference scheme is used to discretize the potential equation

 −ψn+1j+1−2ψn+1j+ψn+1j−1Δx2=qρn+1jϵ0, (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