 # Conservative discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System

We propose a class of conservative discontinuous Galerkin methods for the Vlasov-Poisson system written as a hyperbolic system using Hermite polynomials in the velocity variable. These schemes are designed to be systematically as accurate as one wants with provable conservation of mass and possibly total energy. Such properties in general are hard to achieve within other numerical method frameworks for simulating the Vlasov-Poisson system. The proposed scheme employs discontinuous Galerkin discretization for both the Vlasov and the Poisson equations, resulting in a consistent description of the distribution function and electric field. Numerical simulations are performed to verify the order of accuracy and conservation properties.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

One of the simplest model that is currently applied in plasma physics simulations is the Vlasov-Poisson system. Such a system describes the evolution of charged particles (electrons and ions) under the effects of a self-consistent electric field. For each species labelled , the unknown , depending on the time , the position , and the velocity , represents the distribution function of particles in phase space. This model can be used for the study of beam propagation or a collisionless plasma

 (1.1) ⎧⎪⎨⎪⎩∂fα∂t+v⋅∇xfα+qαmαE⋅∇vfα=0,fα(t=0)=fα,0,

where represents the single charge and represents the mass of one particle , whereas the electric field satisfies the Poisson equation

 (1.2) −4πϵ0ΔxΦ=∑αqαnα,withnα=∫R3fαdv,

where is the vacuum permittivity. On the one hand, for a smooth and nonnegative initial data , the solution to (1.1) remains smooth and nonnegative for all . On the other hand, for any function , we have

 ddt∫Rd×Rdβ(fα(t))dxdv=0,∀t∈R+,

which leads to the conservation of mass, norms, for and kinetic entropy,

 H(t):=∫R6fα(t)ln(f(t))dxdv=H(0),∀t≥0.

We also get the conservation of momentum

 ∑α∫R6mαvfα(t)dxdv=∑α∫R6mαvfα,0dxdv

and total energy

 E(t):=∑αmα2∫R6fα(t)∥v∥2dxdv+2πϵ0∫R3∥E∥2dx=E(0),∀t≥0.

The numerical resolution of the Vlasov-Poisson system (1.1)-(1.2) is often performed by particle methods (PIC) which consist in approximating the plasma by a finite number of macro-particles. The trajectories of these particles are computed from the characteristic curves given by the Vlasov equation, whereas self-consistent fields are computed on a mesh of the physical space. This method allows to obtain satisfying results with a few number of particles (we refer to Birdsall-Langdon for more details ). However, it suffers from poor accuracy since the numerical noise only decreases in when the number of particles

is increased, which is not satisfying for some specific problems. Therefore, different approaches, discretizing the Vlasov equation on a mesh of phase space, have been proposed. Among them, the Fourier-Fourier transform is based on a Fast Fourier Transform of the distribution function in phase space, but this method is only valid for periodic boundary conditions



. Consequently, for non periodic boundary conditions, Gibbs oscillations form at the boundary of the grid and become source of spurious oscillations which propagate into the distribution function. Later semi-Lagrangian methods, which consist in computing the distribution function at each grid point by following the characteristic curves backward, were also used. To compute the origin of the characteristics, a high order interpolation method is needed. E. Sonnendrücker

proposed the cubic spline reconstruction which gives very good results , but the use of spline interpolation destroys the local character of the reconstruction. Nakamura and Yabe also presented the Cubic Interpolated Propagation (CIP) method based on the approximation of the gradients of the distribution function in order to use a Hermite interpolation . This method is very expensive in memory computation since it needs the storage of , , and . Another scheme for the Vlasov equation is the Flux Corrected Transport (FCT)  or more recently finite volume methods [16, 17, 13, 30]: the basic idea is to compute the average of the Vlasov equation solution in each cell of the phase space grid by a conservative method. In the sixties, a variant has been proposed: rather to discretize the function in velocity space, Galerkin methods with a “small” finite set of orthogonal polynomials are used [1, 24]. More recently, in [33, 22] and , the authors have shown the merit to use rescaled orthogonal basis like the so-called scaled Hermite basis. In , the authors show that, according to some symmetry properties of the Hermite weights, the numerical method conserves certain integral quantities like the number of particles, the momentum, the energy or the

-norm. These properties lead to algorithms that can provide long time numerical stability while using a small set of basis functions. The aim of the present work is to apply these techniques in order to handle with the velocity space and design conservative discretization with a reduced degree of freedom.

For the space discretization, we adopt the point of view of discontinuous Galerkin methods. These methods are similar to finite element methods but use discontinuous polynomials and are particularly well suited to handle complicated boundaries which may arise in many practical applications. Moreover, their local construction endows the methods with good local conservation properties without sacrificing the order of accuracy. Furthermore, they are extremely flexible in handling -adaptivity, the boundary conditions are imposed weakly and the mass matrices are block-diagonal which results in very efficient time-stepping algorithms in the context of time-dependent problems, as it is the case here. Among the computational works, we mention the use of discontinuous Galerkin schemes for the Vlasov-Poisson system in [11, 20]; and Vlasov-Maxwell . Several theoretical works have been performed to analyze a family of semi-discrete discontinuous Galerkin schemes for the Vlasov-Poisson system with periodic boundary conditions, for the one  and multi-dimensional cases  , but also for the Vlasov-Maxwell system 

. The authors show optimal error estimates for both the distribution function and the electrostatic field, and they study the conservation properties of the proposed schemes. Let us emphasize that these former works apply a discontinuous Galerkin method using a phase space mesh, whereas here we adopt this approach only in physical space. Furthermore, we propose to modify the fluxes to ensure conservation of mass, momentum and total energy. We prove these conservations for the semi-discrete and fully discrete cases (first and second order time discretizations).

This paper is organized as follows: in the first part (Section 2), we briefly reformulate the Vlasov equation using the Hermite basis in velocity, recalling some properties of the solution. Then, in Section 3, we present a variant of the discontinuous Galerkin method for the space discretization such that mass, momentum and energy are preserved for semi-discrete and fully discrete systems. In the last section (Section 4), we present numerical results for Landau damping, two stream instability and bump on the tail problems to illustrate accuracy and conservations of our discretization technique.

## 2. Hermite spectral form of the Vlasov equation

For simplicity, we now set all the physical constants to one and consider the one dimensional Vlasov-Poisson system for a single species with periodic boundary conditions in space, where the Vlasov equation reads

 (2.1) ∂f∂t+v∂f∂x+E∂f∂v=0

with , position and velocity . Moreover, the self-consistent electric field is determined by the Poisson equation

 (2.2) ∂E∂x=ρ−ρ0,

where the density is given by

 ρ(t,x)=∫Rf(t,x,v)dv,t≥0,x∈(0,L)

and is a constant ensuring the quasi-neutrality condition of the plasma

 (2.3) ∫L0(ρ−ρ0)dx=0.

### 2.1. Hermite basis

To discretize the velocity variable, our approach is based on spectral methods, where we expand the velocity part of the distribution function in basis functions (typically Fourier or Hermite), leading to a truncated set of moment equations for the expansion coefficients. Actually, spectral methods are commonly used to approximate the solution to partial differential equations

[6, 21] and in particular the Vlasov-Poisson system [28, 25, 14, 22, 29].

For instance, in their seminal work , M. Shoucri and G. Knorr applied Chebyshev polynomials, whereas in [25, 14] the authors used Fourier basis but these methods does not conserve neither momentum, nor total energy. Later, J. P. Holloway and J. W. Schumer [22, 29] applied Hermite functions. Indeed, the product of Hermite polynomials and a Gaussian, seems to be a natural choice for Maxwellian-type velocity profiles . Moreover, compared with other classical polynomials, they possess a fairly simple expression for their derivatives and allow to get conservation of mass, momentum and total energy. Recently, these methods generate a new interest leading to new techniques to improve their efficiency [26, 7, 9, 27, 8, 12] .

Let us summarize the main characteristics of this approach. The solution is approximated by a finite sum which corresponds to a truncation of a series

 (2.4) f(t,x,v)=NH−1∑n=0Cn(t,x)Ψn(v),

where is the number of modes. The issue is then to determine a well-suited class of basis functions and to find the expansion coefficients . Here, we have chosen the following basis of normalized scaled asymmetrically weighted Hermite functions:

 (2.5) Ψn(v)=Hn(vvth)e−v2/2v2th√2π,

where corresponds to the scaling velocity and we have set , and for , has the following recursive relation

 (2.6) √nHn(ξ)=ξHn−1(ξ)−√n−1Hn−2(ξ),∀n≥1.

The Hermite basis (2.5) has the following properties

 (2.7) ∫RΨn(ξ)Ψm(ξ)eξ2/2√2πdξ=δn,m,

and is the Kronecker delta function. With those properties, one can substitute the expression (2.4) for into the Vlasov equation (2.1) and using the orthogonality property (2.7), it yields an evolution equation for each mode ,

 (2.8) ∂Cn∂t+vth(√n+1∂Cn+1∂x+√n∂Cn−1∂x)−√nvthECn−1=0.

Meanwhile, we first observe that the density satisfies

 ρ=∫Rfdv=vthC0,

then the Poisson equation becomes

 (2.9) ∂E∂x=vthC0−ρ0.

Observe that if we take in the expression (2.4), we get an infinite system (2.8)-(2.9) of equations for and , which is formally equivalent to the Vlasov-Poisson system (2.1)-(2.2).

### 2.2. Conservation properties

There are different kinds of Hermite approximations based on the choice of Hermite functions and weights. In [22, 29], J. W. Schumer and J. P. Holloway have discussed precisely the different choices of orthogonal Hermite basis functions depending on the form of their weight functions. If the basis is asymmetrically weighted, hence mass, momentum and total energy are preserved but not the norm. In the case of symmetrically weighted basis functions, particles number, total energy (for odd) or momentum (for even), as well as -norm are conserved and the numerical stability is shown to be better. Moreover, it has been shown in  and later in  that it is required to introduce a velocity scaling factor to make the method more accurate and stable. Here, the basis is asymmetrically weighted, hence we prove the following result for the truncated system (2.8)-(2.9).

###### Proposition 2.1.

For any , consider the distribution function given by the truncated series

 f(t,x,v)=NH−1∑n=0Cn(t,x)Ψn(v),

where is a solution to the Vlasov-Poisson system written as (2.8)-(2.9). Then mass, momentum and total energy are conserved, that is,

 vk+1thddt∫L0Ckdx=0,k=0,1

and for the total energy,

 E(t):=12∫L0v3th(√2C2+C0)+|E|2dx=E(0).
###### Proof.

We consider the first three equations of (2.8) given by

 (2.10) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂C0∂t+vth∂C1∂x=0,∂C1∂t+vth(√2∂C2∂x+∂C0∂x)−1vthEC0=0,∂C2∂t+vth(√3∂C3∂x+√2∂C1∂x)−√2vthEC1=0,

which will be related to the conservation of mass, momentum and energy. Thus, let us take a look at the conservation properties from (2.10), by considering periodic or compact support boundary conditions.

First the total mass is defined as

 ∫L0∫∞−∞f(t,x,v)dvdx=vth∫L0C0(t,x)dx,

hence the conservation of mass directly comes from (2.10) by integrating the first equation with respect to .

Then we define the momentum as

 ∫L0∫∞−∞vf(t,x,v)dvdx=v2th∫L0C1(t,x)dx,

due to the second equation (2.10) and using the Poisson equation (2.9), we get

 v2thddt∫L0C1dx=vth∫L0EC0dx=∫L0E(∂E∂x+ρ0)dx=ρ0∫L0Edx=0,

so that the conservation of momentum is obtained.

Finally to prove the conservation of total energy , we compute the variation of the kinetic energy defined as

 12∫L0∫∞−∞v2f(t,x,v)dxdv=v3th2∫L0(√2C2+C0)dx.

Thus, combining the first and third equations in (2.10) and integrating over , we get

 (2.11) v3th2ddt∫L0(√2C2+C0)dx=v2th∫L0EC1dx.

On the other hand, multiplying the first equation in (2.10) by and integrating over , it yields

 vth∫L0∂C0∂tΦdx=−v2th∫L0∂C1∂xΦdx=−v2th∫L0C1Edx.

Using the Poisson equation (2.9), we have , hence

 12ddt∫L0|E|2dx=∫L0∂∂t(−∂2Φ∂x2)Φdx=−v2th∫L0C1Edx.

Combining this latter equality with (2.11), it yields the conservation of total energy

 12ddt∫L0v3th(√2C2+C0)+|E|2dx=0.

###### Remark 2.2.

The natural space associated to our asymmetrically weighted basis is

 L2ω={u:R↦R:∫R|u(v)|2e|v|2/2v2thdv<∞}.

Unfortunately there is no estimate of the associated norm for the solution to the Vlasov-Poisson system (2.1)-(2.2), hence there is no warranty to get such estimate for the system (2.8)-(2.9). It is worth to mention that the symmetrically weighted basis with the norm would be a good choice, but it would affect the efficiency of the conservation algorithm which will proposed in the Section 3 (see Remark 3.2).

### 2.3. Filtering technique

Filtering is a common procedure to reduce the effects of the Gibbs phenomenon inherent to spectral methods . The filter will consist in multiplying some spectral coefficients in (2.4) by a scaling factor in order to reduce the amplitude of high frequencies, for any ,

 ˜Cn=Cnσ(nNH).

Here, we simply apply a filter, called Hou-Li’s filter, already proposed in  for Fourier spectral method, which reads

 σ(s)={1,if\,\,0≤|s|≤2/3,exp(−β|s|β),if \,\,|s|>2/3,

where the coefficient is chosen as . Actually this filter achieved great success in Fourier  and Hermite  spectral methods. We again refer to  for a detailled discussion on the application to the Vlasov-Poisson system using a Hermite spectral method.

###### Remark 2.3.

Observe that the filter is applied only when , hence the filtering process does not modify the coefficients , that is, conservations of mass, momentum and total energy are not affected.

## 3. Discontinuous Galerkin scheme

This section is devoted to the space and time discretizations of the Vlasov-Poisson system written in the form (2.8)-(2.9). We apply a local discontinuous Galerkin method for space and introduce a slight modification of the fluxes to preserve mass, momentum and local energy for this semi-discrete system. Then, we focus on the time discretization and prove the conservation for the fully-discrete case.

### 3.1. Semi-discrete discontinuous Galerkin scheme

In this section, we will define the discontinuous Galerkin (DG) scheme for the Vlasov-Poisson system with Hermite spectral basis in velocity (2.8)-(2.9).

Let us first introduce some notations and start with , a partition of . Here , . Each element is denoted as with its length , and .

Given any non-negative integer , we define a finite dimensional discrete piece-wise polynomial space

where the local space consists of polynomials of degree at most on the interval . We further denote the jump and the average of at defined as,

 [u]i+12=u(x+i+12)−u(x−i+12)and{u}i+12=12(u(x+i+12)+u(x−i+12)),∀i,

where . We also denote

 ⎧⎪⎨⎪⎩ui+12=u(xi+12),u±i+12=u(x±i+12).

From these notations, we apply a semi-discrete discontinuous Galerkin method for (2.8) as follows. We look for an approximation , such that for any , we have

 (3.1) ddt∫IjCn,hϕndx+ajn(gn,ϕn)−√nvth∫IjECn−1,hϕndx=0,0≤n≤NH−1,

where is defined as

 (3.2) ⎧⎪ ⎪⎨⎪ ⎪⎩ajn(gn,ϕn)=−∫Ijgnϕ′ndx+^gn,j+12ϕ−n,j+12−^gn,j−12ϕ+n,j−12,gn=vth(√n+1Cn+1,h+√nCn−1,h).

The numerical flux in (3.2) could be taken as any consistent numerical fluxes. As an example, the global Lax-Friedrichs flux is used, which is defined as

 (3.3) ^gn=12[g−n+g+n−α(C+n,h−C−n,h)],

with the numerical viscosity coefficient .

For the Poisson equation (2.9), we introduce the potential function , such that

Hence we get the one dimensional Poisson equation

 −∂2Φ∂x2=vthC0−ρ0,

for which we also define a local discontinuous Galerkin scheme: we look for a couple , such that for any and , we have

 (3.4) ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩+∫IjΦhη′dx−^Φh,j+12η−j+12+^Φh,j−12η+j−12=∫Ehηdx,−∫IjEhζ′dx+^Eh,j+12ζ−j+12−^Eh,j−12ζ+j−12=∫(vthC0,h−ρ0)ζdx,

where the numerical fluxes and in (3.4) here are taken as in [3, 2],

 (3.5) {^Φh={Φh},^Eh={Eh}−β[Φh],

with either being a positive constant or a constant multiplying by and is the mesh size.

### 3.2. Conservation properties for the semi-discrete scheme

The local discontinuous Galerkin method presented in the previous section naturally preserves the mass since the equation on only contains a convective term. Unfortunately, the momentum and total energy are not conserved due to the contribution of the source terms in (3.1). Therefore, we propose a slight modification of the fluxes for the unknowns in order to ensure the right conservations.

In this subsection, without any confusion, we will drop the subindex for simplicity. The aim is to modify the scheme (3.1)-(3.5) in order to ensure the conservation of mass, momentum and energy without deteriorating the order of accuracy. Thus, we follow the proof of Proposition 2.1 and consider the scheme (3.1) for combined with (3.4), by adding two residual terms to ensure the conservation of momentum and total energy

 (3.6) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩ddt∫IjC1ϕ1dx+aj1(g1,ϕ1)−1vth∫IjEC0ϕ1dx=⟨rj1,ϕ1⟩,ddt∫IjC2ϕ2dx+aj2(g2,ϕ2)−√2vth∫IjEC1ϕ2dx=⟨rj2,ϕ2⟩,

and the residual terms are given by

 (3.7) ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩⟨rj1,ϕ1⟩=−β2v2th([Φ][E]j−12ϕ+1,j−12+[Φ][E]j+12ϕ−1,j+12),⟨rj2,ϕ2⟩=1√2vth(({C1}−^C1)[Φ]j−12ϕ+2,j−12+({C1}−^C1)[Φ]j+12ϕ−2,j+12).

Apparently, these residues are given as products of two cell interface jump terms, which both are of high order, hence it should not deteriorate the order of accuracy. Therefore, we prove that this choice will ensure the discrete conservation properties of the scheme.

###### Proposition 3.1.

For any , we consider the solution to the system (3.1)-(3.5) combined with the residual terms (3.6)-(3.7). Then, discrete mass and momentum are conserved

 vk+1thddt∫L0Ckdx=0,k=0,1

and the discrete total energy, defined as

 Eh(t):=12∫L0v3th(√2C2+C0)+|E|2dx+12β∑j[Φ]2j−12,

is also preserved with respect to time.

###### Proof.

On the one hand, the conservation of mass easily follows by choosing in (3.1) for and summing over , hence since there is no source term we get

 ∑jddt∫IjC0dx=0.

On the other hand, for the conservation of momentum, taking in the modified scheme (3.6) and summing over , we have

 (3.8) ddt∫L0C1dx=1vth∫L0EC0dx+∑j⟨rj1,1⟩,

where the contribution of the residual term, using periodic boundary conditions, gives

 (3.9) ∑j⟨rj1,1⟩=−∑jβ2v2th([Φ][E]j−12+[Φ][E]j+12)=−βv2th∑j[Φ][E]j−12.

Then we compute the contribution of the first source term on the right hand side of (3.8). Choosing in the first equation of (3.4), we get that

 ∫IjEdx=−^Φj+12+^Φj−12,

whereas taking in the second equation of (3.4), we have

 −∫IjE∂E∂xdx+^Ej+12E−j+12−^Ej−12E+j−12=∫Ij(vthC0−ρ0)Edx,

hence adding the latter both equalities, it yields

 vth∫IjEC0dx =−12∫Ij∂E2∂xdx+^Ej+12E−j+12−^Ej−12E+j−12+ρ0∫IjEdx =−12(E2)−j+12+12(E2)+j−12+^Ej+12E−j+12−^Ej−12E+j−12+ρ0(^Φj+12−^Φj−12).

Summing on , we obtain

 (3.10) vth∫L0EC0dx=∑j(12[E2]−^E[E])j−12=∑j({E}−^E)j−12[E]j−12.

Finally, using the definition of the flux in (3.5), we have

 ({E}−^E)j−12[E]j−12=β[Φ][E]j−12,

hence gathering (3.9) and (3.10) in (3.8), we get the conservation of momentum

 ddt∫L0C1dx=0.

Now let us show the conservation of total energy. Choosing in the modified scheme (3.6) and summing over , we get

 (3.11) ddt∫L0C2dx=√2vth∫L0EC1dx+∑j⟨rj2,1⟩,

where the contribution of the residual term, using periodic boundary conditions, gives

 (3.12) ∑j⟨rj2,1⟩=√2vth∑j(({C1}−^C1)[Φ])j−12.

Let us evaluate the first term on the right hand side of (3.11). We consider the first equation of (3.4) and take ,

 ∫IjEC1dx=∫IjΦ∂C1∂xdx−^Φj+12(C1)−j+12+^Φj−12(C1)+j−12,

whereas in (3.1), we choose ,

 1vth∫Ij∂C0∂tΦdx=∫IjC1∂Φ∂xdx−^C1,j+12Φ−j+12+^C1,j−12Φ+j−12.

Adding the latter two equalities, it yields

 ∫IjEC1+1vth∂C0∂tΦdx =∫Ij∂(ΦC1)∂xdx−^Φj+12(C1)−j+12+^Φj−12(C1)+j−12−^C1,j+12Φ−j+12+^C1,j−12Φ+j−12 =(ΦC1)−j+12−(ΦC1)+j−12−^Φj+12(C1)−j+12+^Φj−12(C1)+j−12−^C1,j+12Φ−j+12+^C1,j−12Φ+j−12.

Summing on and using periodic boundary conditions, we get

 ∫L0[EC1+1vth∂C0∂tΦ]dx=−∑j([ΦC1]−^Φ[C1]−^C1[Φ])j−12.

To evaluate the right hand side of the latter term, we use that

 [ΦC1]−^Φ[C1]−^C1[Φ]=({Φ}−^Φ)[C1]+({C1}−^C1)[Φ],

and from the definition of the flux for the Poisson equation (3.5), we get

 (3.13) √2vth∫L0[EC1+1vth∂C0∂tΦ]dx=−√2vth∑j(({C1}−^C1)[Φ])j−12.

Now, it remains to evaluate the second term on the left hand side in (3.13), hence we compute the time derivative of the second equation in (3.4) and choose the test function , it gives

 √2v2th∫Ij∂C0∂tΦdx=√2v3th⎡⎢⎣−∫Ij∂E∂t∂Φ∂xdx+(ˆ∂E∂t)j+12Φ−j+12−(ˆ∂E∂t)j−12Φ+j−12⎤⎥⎦.

Then in the first equation of (3.4), we take and obtain

 −1√2v3thddt∫IjE2dx = −√2