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
where represents the single charge and represents the mass of one particle , whereas the electric field 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
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ückerproposed 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
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
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
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:
where corresponds to the scaling velocity and we have set , and for , has the following recursive relation
The Hermite basis (2.5) has the following properties
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 ,
Meanwhile, we first observe that the density satisfies
then the Poisson equation becomes
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).
We consider the first three equations of (2.8) given by
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
hence the conservation of mass directly comes from (2.10) by integrating the first equation with respect to .
Then we define the momentum as
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
Thus, combining the first and third equations in (2.10) and integrating over , we get
On the other hand, multiplying the first equation in (2.10) by and integrating over , it yields
Using the Poisson equation (2.9), we have , hence
Combining this latter equality with (2.11), it yields the conservation of total energy
The natural space associated to our asymmetrically weighted basis is
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 ,
Here, we simply apply a filter, called Hou-Li’s filter, already proposed in  for Fourier spectral method, which reads
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.
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
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,
where . We also denote
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
where is defined as
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
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
for which we also define a local discontinuous Galerkin scheme: we look for a couple , such that for any and , we have
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
and the residual terms are given by
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.
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
On the other hand, for the conservation of momentum, taking in the modified scheme (3.6) and summing over , we have
where the contribution of the residual term, using periodic boundary conditions, gives
whereas taking in the second equation of (3.4), we have
hence adding the latter both equalities, it yields
Summing on , we obtain
Finally, using the definition of the flux in (3.5), we have
Now let us show the conservation of total energy. Choosing in the modified scheme (3.6) and summing over , we get
where the contribution of the residual term, using periodic boundary conditions, gives
whereas in (3.1), we choose ,
Adding the latter two equalities, it yields
Summing on and using periodic boundary conditions, we get
To evaluate the right hand side of the latter term, we use that
and from the definition of the flux for the Poisson equation (3.5), we get
Then in the first equation of (3.4), we take and obtain