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

06/14/2021 ∙ by Marianne Bessemoulin-Chatard, et al. ∙ 0

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.



There are no comments yet.


page 23

page 25

page 28

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 


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


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

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

We also get the conservation of momentum

and total energy

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


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:

with the corresponding norm, that is

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

Then, since

we obtain

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

We notice that if is such that


the first and third terms cancel each other, yielding

Finally, applying Grönwall’s Lemma gives

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


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 :

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


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


where the density is given by

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

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


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:


where is the positive nonincreasing function given by


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

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


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

the transport term is

and finally the nonlinear term is

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


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

and then the Poisson equation becomes


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 ,


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


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

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,

and for the total energy,


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


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

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

We have

which taking in the second equation of (2.11) yields

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

which gives the conservation of momentum.

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

We have

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


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

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

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

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

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

2.3. Weighted stability of

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

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

Proposition 2.2.

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

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


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

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


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

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

On the other hand, we notice that

Gathering these two identities in (2.14), we get

Applying Young inequality in the second sum, this provides:

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

which yields

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 ,


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

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

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

and then one obtains

Hence, we need to control