# On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system

We study a class of spatial discretizations for the Vlasov-Poisson system written as an hyperbolic system using Hermite polynomials. In particular, we focus on spectral methods and discontinuous Galerkin approximations. To obtain L 2 stability properties, we introduce a new L 2 weighted space, with a time dependent weight. For the Hermite spectral form of the Vlasov-Poisson system, we prove conservation of mass, momentum and total energy, as well as global stability for the weighted L 2 norm. These properties are then discussed for several spatial discretizations. Finally, numerical simulations are performed with the proposed DG/Hermite spectral method to highlight its stability and conservation features.

## Authors

• 1 publication
• 6 publications
• ### Conservative discontinuous Galerkin/Hermite Spectral Method for the Vlasov-Poisson System

We propose a class of conservative discontinuous Galerkin methods for th...
04/06/2020 ∙ by Francis Filbet, et al. ∙ 0

• ### Stability and conservation properties of Hermite-based approximations of the Vlasov-Poisson system

Spectral approximation based on Hermite-Fourier expansion of the Vlasov-...
03/01/2021 ∙ by Daniele Funaro, et al. ∙ 0

• ### Analysis of an exactly mass conserving space-time hybridized discontinuous Galerkin method for the time-dependent Navier–Stokes equations

We introduce and analyze a space-time hybridized discontinuous Galerkin ...
03/24/2021 ∙ by Keegan L. A. Kirk, et al. ∙ 0

• ### Stability of Discontinuous Galerkin Spectral Element Schemes for Wave Propagation when the Coefficient Matrices have Jumps

We use the behavior of the L_2 norm of the solutions of linear hyperboli...
11/23/2020 ∙ by David A Kopriva, et al. ∙ 0

• ### Some continuous and discontinuous Galerkin methods and structure preservation for incompressible flows

In this paper, we present consistent and inconsistent discontinuous Gale...
02/18/2021 ∙ by Xi Chen, et al. ∙ 0

• ### Petrov-Galerkin flux upwinding for mixed mimetic spectral elements, and its application to vorticity stabilisation in shallow water

Upwinded mass fluxes are described and analysed for advection operators ...
04/28/2020 ∙ by David Lee, et al. ∙ 0

• ### Stability and Convergence of Spectral Mixed Discontinuous Galerkin Methods for 3D Linear Elasticity on Anisotropic Geometric Meshes

We consider spectral mixed discontinuous Galerkin finite element discret...
08/13/2019 ∙ by Thomas P. Wihler, et al. ∙ 0

##### 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. This system describes the evolution of charged particles in the case where the only interaction considered is the mean-field force created through electrostatic effects. The system consists in Vlasov equations for phase space density of each particle species of charge and mass

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

coupled to its self-consistent electric field which satisfies the Poisson equation

 (1.2) −4πϵ0ΔxΦ=∑βqβnβ,withnβ=∫Rdfβ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∫R2dG(fβ(t))dxdv=0,∀t∈R+,

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

 H(t):=∫R2dfβ(t)ln(f(t))dxdv=H(0),∀t≥0.

We also get the conservation of momentum

and total energy

 E(t):=∑βmβ2∫R2dfβ(t)∥v∥2dxdv+2πϵ0∫Rd∥E∥2dx=E(0),∀t≥0.

Numerical approximation of the Vlasov-Poisson system has been addressed since the sixties. Particle methods (PIC), consisting in approximating the plasma by a finite number of macro particles, have been widely used [4]. They allow to obtain satisfying results with a few number of particles, but a well-known drawback of this class of methods is their inherent numerical noise which only decreases in when the number of particles increases, preventing from getting an accurate description of the distribution function for some specific applications. To overcome this difficulty, Eulerian solvers, that is methods discretizing the Vlasov equation on a mesh of the phase space, can be considered. Many authors explored their design, and an overview of many different techniques with their pros and cons can be found in [11]. Among them, we can mention finite volume methods [10]

which are a simple and inexpensive option, but in general low order. Fourier-Fourier transform schemes

[19] are based on a Fast Fourier Transform of the distribution function in phase space, but suffer from Gibbs phenomena if other than periodic conditions are considered. Standard finite element methods [31, 32] have also been applied, but may present numerical oscillations when approximating the Vlasov equation. Later, semi-Lagrangian schemes have also been proposed [29]

, consisting in computing the distribution function at each grid point by following the characteristic curves backward. Despite these schemes can achieve high order allowing also for large time steps, they require high order interpolation to compute the origin of the characteristics, destroying the local character of the reconstruction.

In the present article, using Hermite polynomials in the velocity variable, we write the Vlasov-Poisson system (1.1)–(1.2) as an hyperbolic system. This idea of using Galerkin methods with a small finite set of orthogonal polynomials rather than discretizing the distribution function in velocity space goes back to the 60’s [1, 18]. More recently, the merit to use rescaled orthogonal basis like the so-called scaled Hermite basis has been shown [8, 16, 28, 26, 30]. In [16], Holloway formalized two possible approaches. The first one, called symmetrically-weighted (SW), is based on standard Hermite functions as the basis in velocity and as test functions in the Galerkin method. It appears that this SW method cannot simultaneously conserve mass and momentum. It makes up for this deficiency by correctly conserving the norm of the distribution function, ensuring the stability of the method. In the second approach, called asymmetrically-weighted (AW), another set of test functions is used, leading to the simultaneous conservation of mass, momentum and total energy. However, the AW Hermite method does not conserve the norm of the distribution function and is then not numerically stable. In addition, the asymmetric Hermite method exactly solves the spatially uniform problem without any truncation error, a property not shared by either the traditional symmetric Hermite expansion or by finite difference methods. The aim of this work is to present a class of numerical schemes based on the AW Hermite methods and to provide a stability analysis.

In what follows, we consider two types of spatial discretizations for the Vlasov-Poisson system (1.1)-(1.2), written as an hyperbolic system using Hermite polynomials in the velocity variable: spectral methods and a discontinuous Galerkin (DG) method. Concerning spectral methods, the Fourier basis is the natural choice for the spatial discretization when considering periodic boundary conditions. Spectral Galerkin and spectral collocation methods for the AW Fourier-Hermite discretization have been proposed in [8, 21, 24]. In [5], authors study a time implicit AW Fourier-Hermite method allowing exact conservation of charge, momentum and energy, and highlight that for some test cases, this scheme can be significantly more accurate than the PIC method.

For the SW Fourier-Hermite method, a convergence theory was proposed in [25]. In [20], authors study conservation and stability properties of a generalized Hermite-Fourier semi-discretization, including as special cases the SW and AW approaches. Concerning discontinuous Galerkin methods, they are similar to finite elements methods but use discontinuous polynomials and are particularly well-adapted to handling complicated boundaries which may arise in many realistic applications. Due to their local construction, this type of methods provides good local conservation properties without sacrificing the order of accuracy. They were already used for the Vlasov-Poisson system in [14, 6]

. Optimal error estimates and study of the conservation properties of a family of semi-discrete DG schemes for the Vlasov-Poisson system with periodic boundary conditions have been proved for the one

[2] and multi-dimensional [3] cases. In all these works, the DG method is employed using a phase space mesh.

Here, we adopt this approach only in physical space, as in [13], with a Hermite approximation in the velocity variable. In [13], such schemes with discontinuous Galerkin spatial discretization are designed in such a way that conservation of mass, momentum and total energy is rigorously provable.

As mentionned before, the main difficulty is to study the stability of approximations based on asymmetrically-weighted Hermite basis. Indeed, this choice fails to preserve norm of the approximate solution, and therefore to ensure long-time stability of the method. In this framework, the natural space to be considered is

and there is no estimate of the associated norm for the solution to the Vlasov-Poisson system (1.1)-(1.2). To overcome this difficulty, we introduce a new weighted space, with a time-dependent weight, allowing to prove global stability of the solution in this space. Actually, this idea has been already employed in [22, 23] to stabilize Hermite spectral methods for linear diffusion equations and nonlinear convection-diffusion equations in unbounded domains, yielding stability and spectral convergence of the considered methods. Here, we define the weight as

 (1.3) ω(t,v):=(2π)d/2e(α(t)∥v∥)2/2,

where is a nonincreasing positive function which will be designed in such a way that a global stability estimate can be established in the following weighted space:

 L2ω(t):={u:R2d→R:∫R2d|u(x,v)|2ω(t,v)dxdv<+∞},

with the corresponding norm, that is

 ∥u∥2ω(t)=(2π)d/2∫R2d|u(x,v)|2e(α(t)∥v∥)2/2dxdv.

Let us now determine the function . To do so, we compute the time derivative of , being the solution of (1.1). Using the Vlasov equation (1.1) and the definition (1.3) of the weight , one has

 12ddt∥fβ(t)∥2ω(t)= −∫R2dfβ(v⋅∇xfβ+qβmβE⋅∇vfβ)ωdxdv +12∫R2dαα′∥v∥2f2βωdvdx.

Then, since

 ∫RdfβE⋅∇vfβωdv=−12∫Rdα2f2βE⋅vωdv,

we obtain

 12ddt∥fβ(t)∥2ω(t)=12∫R2df2β(qβmβα2E⋅v+αα′∥v∥2)ωdxdv.

Applying now Young inequality on the first term, we get for ,

 12ddt∥fβ(t)∥2ω(t)≤12∫R2df2β(γ2q2βm2βα4∥E∥2∞∥v∥2+12γ+αα′∥v∥2)ωdxdv.

We notice that if is such that

 (1.4)

the first and third terms cancel each other, yielding

 12ddt∥fβ(t)∥2ω(t)≤14γ∥fβ(t)∥2ω(t).

Finally, applying Grönwall’s Lemma gives

 ∥fβ(t)∥ω(t)≤∥fβ,0∥ωet/4γ.

It now remains to define satisfying (1.4). We remark that

 (1.5) α(t):=α0(1+γq2βm2β∫t0∥E(s)∥2∞ds)−1/2

is a suitable choice, where is the initial value of at and corresponds to the initial scale of the distribution function, whereas is a free parameter. This function is positive, nonincreasing and the parameter will practically be chosen sufficiently small to ensure that does not decrease too fast towards as goes to infinity.

Summarizing, we have established the following result.

###### Proposition 1.1.

Assuming that the initial data belongs to , then the solution to (1.1) satisfies for all :

 ∥fβ(t)∥ω(t)≤∥fβ,0∥ω0et/4γ,

where is the constant appearing in the definition (1.5) of .

Notice that the weight depends on the solution to the Vlasov-Poisson, hence it is mandatory to control the norm of in order that this estimate makes sense. In the next section, we introduce the formulation of the Vlasov equation using the Hermite basis in velocity, and prove the analogous of Proposition 1.1 for the obtained system. Then in Section 3, we discuss conservation and stability properties for a class of spatial discretizations including Fourier spectral method and discontinuous Galerkin approximations. In Section 4, we introduce the time discretization that will be used to compute numerical results with the discontinuous Galerkin method. Finally in Section 5 we present numerical results for two stream instability, bump-on-tail problem and ion acoustic wave to highlight conservations and stability of our discretization.

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

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

### 2.1. Hermite basis

We approximate the solution of (2.1)–(2.2) by a finite sum which corresponds to a truncation of a series

 (2.3) fNH(t,x,v)=NH−1∑n=0Cn(t,x)Ψn(t,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, our aim is to obtain a control on a norm of , in the spirit of that established in Proposition 1.1. We then choose the following basis of normalized scaled time dependent asymmetrically weighted Hermite functions:

 (2.4) Ψn(t,v)=α(t)Hn(α(t)v)e−(α(t)v)2/2√2π,

where is the positive nonincreasing function given by

 (2.5) α(t)=α0(1+γ∫t0∥ENH(s)∥2∞ds)−1/2,

with an approximation of the electric field that will be defined below. Functions are the Hermite polynomials defined by , and for , has the following recursive relation

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

Let us also emphasize that for all , and that the Hermite basis (2.4) satisfies the following orthogonality property

 (2.6) ∫RΨn(t,v)Hm(α(t)v)dv=δn,m,

where is the Kronecker delta function.

Now the objective is to obtain an evolution equation for each mode by substituting the expression (2.3) for into the Vlasov equation (2.1) and using the orthogonality property (2.6). Thanks to the definition (2.4) of and the properties of , we compute the different terms of (2.1). The time derivative term is given by

 ∂tfNH=NH−1∑n=0(∂tCnΨn−Cnα′α(nΨn+√(n+1)(n+2)Ψn+2)),

the transport term is

 v∂xfNH=NH−1∑n=0∂xCn1α(√n+1Ψn+1+√nΨn−1),

and finally the nonlinear term is

 ENH∂vfNH=−NH−1∑n=0ENHCnα√n+1Ψn+1.

Then taking as test function and using orthogonality property (2.6), we obtain an evolution equation for each mode , :

 (2.7) ∂tCn−α′α(nCn+√(n−1)nCn−2)+1α(√n∂xCn−1+√n+1∂xCn+1)−ENHα√nCn−1=0,

with the understanding that for and . Meanwhile, we first observe that the density satisfies

 ρNH=∫RfNHdv=C0,

and then the Poisson equation becomes

 (2.8) ∂ENH∂x=C0−ρ0.

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

In what follows, we rather consider a generic weak formulation of (2.7)–(2.8). Indeed, this will allow us to straightforwardly apply the obtained results to spatial discretizations with spectral methods (see Section 3.1).

Let be a subspace of . We look for such that and for all ,

 ddt∫L0Cnφndx −α′α∫L0(nCn+√(n−1)nCn−2)φndx (2.9) −1α∫L0(√nCn−1+√n+1Cn+1)φ′ndx −α√n∫IjENHCn−1φndx=0,0≤n≤NH−1,

and for a couple such that , and for all and , we have

 (2.10) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∫L0ΦNHη′dx=∫L0ENHηdx,−∫L0ENHζ′dx=∫L0(C0−ρ0)ζdx.

In the rest of this section, we consider , and then (2.9)–(2.10) is the weak formulation of (2.7)–(2.8).

### 2.2. Conservation properties

###### Proposition 2.1.

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

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

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

 ddt∫L0Ckαkdx=0,k=0,1,

and for the total energy,

 ENH(t):=12∫L0(1α2(√2C2+C0)+|ENH|2)dx=ENH(0).
###### Proof.

We consider the first three equations of (2.9): for all ,

 (2.11) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ddt∫L0C0φ0dx−1α∫L0C1φ′0dx=0,ddt∫L0C1φ1dx−α′α∫L0C1φ1dx−1α∫L0(C0+√2C2)φ′1dx−α∫L0ENHC0φ1dx=0,ddt∫L0C2φ2dx−α′α∫L0(2C2+√2C0)φ2dx−1α∫L0(√2C1+√3C3)φ′2dx−√2α∫L0ENHC1φ2dx=0,

which will be related to the conservation of mass, momentum and energy. Thus, let us look at the conservation properties from (2.11).

First the total mass is defined as

 ∫L0∫RfNH(t,x,v)dvdx=∫L0C0(t,x)dx,

hence the conservation of mass directly comes from (2.11) by taking as test function in the first equation.

Then the momentum is given by

 ∫L0∫RvfNH(t,x,v)dvdx=∫L0C1(t,x)α(t)dx.

We have

 ddt∫L0C1αdx=∫L0(∂tC1α−α′α2C1)dx,

which taking in the second equation of (2.11) yields

 ddt∫L0C1αdx=∫L0ENHC0dx.

Now, choosing in the second equation of (2.10) gives

 ddt∫L0C1αdx=ρ0∫L0ENHdx−∫L0ENH∂xENHdx=ρ0∫L0ENHdx=0,

which gives the conservation of momentum.

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

 EKNH(t)=12∫L0∫Rv2fNH(t,x,v)dxdv=121α2∫L0(√2C2+C0)dx.

We have

 dEKNH(t)dt=12α2∫L0(√2∂tC2+∂tC0)dx−α′α3∫L0(√2C2+C0)dx.

Thus, combining the first equation of (2.9) with and the third equation of (2.9) with , we get

 (2.12) dEKNH(t)dt=1α∫L0ENHC1dx.

Finally, taking in the first equation of (2.10), we obtain

 dEKNH(t)dt=1α∫L0ΦNH∂xC1dx.

We now compute the time derivative of the potential energy defined by

 EPNH(t):=12∫L0|ENH|2dx.

Using as test function in the first equation of (2.10) and an integration by parts, we have

 dEPNH(t)dt=∫L0∂tENHENHdx=∫L0ΦNH∂x∂tENHdx=−∫L0∂xΦNH∂tENHdx.

Now, choosing in the time derivative of the second equation of (2.10), and then in the first equation of (2.9), we finally get

 dEPNH(t)dt=∫L0∂tC0ΦNHdx=1α∫L0C1∂xΦNHdx=−1α∫L0∂xC1ΦNHdx=−dEKNH(t)dt.

This concludes the proof since the total energy is the sum of the kinetic and potential energies.

### 2.3. Weighted L2 stability of fNH

We now establish the discrete counterpart of Proposition 1.1, namely the stability in an appropriately weighted norm. More precisely, with defined by (2.5), we consider the weight , and denote by the corresponding weighted norm. Using the definition (2.3) of , we have

 ∥fNH(t)∥2ω(t)=NH−1∑n,m=0∫L0∫RαCnCmΨm(t,v)Hn(α(t)v)dvdx,

which gives, using orthogonality property (2.6), the following expression for the weighted norm of :

 (2.13) ∥fNH(t)∥2ω(t)=NH−1∑n=0α∫L0C2ndx.
###### Proposition 2.2.

Assuming that , the distribution function given by the truncated series (2.3) satisfies the following stability estimate:

 ∥fNH(t)∥ω(t)≤∥fNH(0)∥ω(0)et/4γ,∀t≥0,

where is the constant appearing in the definition (2.5) of .

###### Proof.

We compute the time derivative of defined by (2.13):

 ddt∥fNH(t)∥2ω(t)=NH−1∑n=0(α∫L0∂tCnCndx+α′2∫L0C2ndx).

Thanks to equation (2.9) with , we then have to estimate

 12ddtNH−1∑n=0α(t)∫L0Cn(t,x)2dx= NH−1∑n=0∫L0α′[(n+12)C2n+√n(n−1)CnCn−2]dx +NH−1∑n=0∫L0(√nCn−1+√n+1Cn+1)∂xCndx (2.14) +NH−1∑n=0∫L0α2√nENHCnCn−1dx.

First of all, the transport term vanishes. Indeed, reindexing the sum over and using that , we have

 NH−1∑n=0∫L0(√nCn−1+√n+1Cn+1)∂xCndx=NH−1∑n=0∫L0√n∂x(Cn−1Cn)dx=0.

Then on the one hand, using again that for all and , we have

 NH−1∑n=0[(n+12)C2n+√n(n−1)CnCn−2]=12NH+1∑n=1(√nCn+√n−1Cn−2)2.

On the other hand, we notice that

 NH−1∑n=0√nCnCn−1=12NH+1∑n=1Cn−1(√nCn+√n−2Cn−2).

Gathering these two identities in (2.14), we get

 12ddtNH−1∑n=0α∫L0C2ndx= +12NH+1∑n=1∫L0ENHα2Cn−1(√nCn+√n−1Cn−2)dx.

Applying Young inequality in the second sum, this provides:

 12ddtNH−1∑n=0α∫L0C2ndx≤ +12NH+1∑n=1∫L0(γ2∥ENH∥2∞α3(√nCn+√n−1Cn−2)2+12γαC2n−1)dx.

Yet, using definition (2.5) of , we have that

 α′(t)=−γ2∥ENH(t)∥2∞α(t)3,

which yields

 12ddtNH−1∑n=0α∫L0C2ndx≤ 14γNH+1∑n=1α∫L0C2n−1dx ≤ 14γNH−1∑n=0α∫L0C2ndx.

Using Grönwall’s lemma, this gives the expected result. ∎

### 2.4. Control of α

Using the definition (2.13) of , the stability result established in Proposition 2.2 gives a stability result in for the coefficients , provided that is bounded from below by a positive constant. Due to definition (2.5) of , this is achieved as soon as is controlled. This result is given in the following proposition.

###### Proposition 2.3.

Under assumptions of Proposition 2.2, the solution of (2.10) satisfies that for all , there exists a constant depending on such that for all ,

 ∥ENH(t)∥∞≤CT.
###### Proof.

Since we are in one space dimension, by Sobolev inequality, there exists a constant such that for all ,

 ∥ENH(t)∥2∞≤c∥ENH(t)∥2H1.

Moreover, since , Poincaré-Wirtinger inequality applies and then there exists such that

 ∥ENH(t)∥2∞≤c′∥∂xENH(t)∥22.

Taking in the second equation of (2.10) and applying Cauchy-Schwarz inequality gives

 ∥∂xENH(t)∥22=∫L0∂xENH(C0−ρ0)dx=∫L0∂xENHC0dx≤∥∂xENH(t)∥2∥C0(t)∥2,

and then one obtains

 ∥ENH(t)∥2∞≤c′∥C0(t)∥22.

Hence, we need to control