Study of conservation and recurrence of RungeKutta discontinuous Galerkin schemes for VlasovPoisson systems
Abstract
In this paper we consider RungeKutta discontinuous Galerkin (RKDG) schemes for VlasovPoisson systems that model collisionless plasmas. Onedimensional systems are emphasized. The RKDG method, originally devised to solve conservation laws, is seen to have excellent conservation properties, be readily designed for arbitrary order of accuracy, and capable of being used with a positivitypreserving limiter that guarantees positivity of the distribution functions. The RKDG solver for the Vlasov equation is the main focus, while the electric field is obtained through the classical representation by Green’s function for the Poisson equation. A rigorous study of recurrence of the DG methods is presented by Fourier analysis, and the impact of different polynomial spaces and the positivitypreserving limiters on the quality of the solutions is ascertained. Several benchmark test problems, such as Landau damping, the twostream instability, and the KEEN (Kinetic Electro static Electron Nonlinear) wave, are given.
Keywords: VlasovPoisson, discontinuous Galerkin methods, recurrence, positivitypreserving, BGK mode, KEEN wave.
1 Introduction
The VlasovPoisson (VP) system is an important equation for modeling collisionless plasmas, one that possesses computational difficulties of more complete kinetic theories. Thus, it serves as an important test bed for algorithm development. The VP system describes the evolution of , the probability distribution function () for finding an electron (at position with velocity at time ) with a uniform background of fixed ions under a selfconsistent electrostatic field. In particular, the nondimensionalized VP system (with time scaled by the inverse plasma frequency and length scaled by the Debye length ) is given by
(1)  
Here the domain , where can be either a finite domain or . The boundary conditions for the above systems are summarized as follows: as or . If is finite, then we can impose either inflow boundary conditions with on , where is the outward normal vector, or more simply impose periodic boundary conditions. For simplicity of discussion, in this paper, we will always assume periodicity in . Also, we add that when the VP system is applied to plasmas the total charge neutrality condition, , is imposed.
The following physical quantities associated with this system are related to its conservation properties:
charge density  
momentum density  (2)  
kinetic energy density  
electrostatic energy density 
Indeed, it is wellknown that the VP system conserves the total electron charge , momentum , and energy . Moreover, any functional of the form is a constant of motion. In particular, this includes the th order invariant and the entropy . Sometimes the functional is also called the enstrophy, and all of these invariants are called Casimir invariants (see, e.g., [34]).
The VP system has been studied extensively for the simulation of collisionless plasmas. Popular numerical approaches include ParticleInCell (PIC) methods [6, 24], Lagrangian particle methods [4, 17], semiLagrangian methods [8, 45], the WENO method coupled with Fourier collocation [58], finite volume (flux balance) methods [7, 18, 19], FourierFourier spectral methods [27, 28], continuous finite element methods [49, 50], among many others. In this paper, we will focus on the discontinuous Galerkin (DG) method to solve the VP system. The original DG method was introduced by Reed and Hill [42] for neutron transport. Lesaint and Raviart [32] performed the first convergence study for the original DG method. Cockburn and Shu in a series of papers [14, 13, 12, 11, 15] developed the RungeKutta DG (RKDG) method for hyperbolic equations. The RKDG methods have been used to simulate the VP system in plasmas by Heath, Gamba, Morrison and Michler [23, 22] and for the gravitational infinite homogeneous stellar system by Cheng and Gamba [9]. Theoretical aspects about stability, accuracy and conservation of those methods are discussed in [22, 23] and more recently in [3, 2] for energy conserving schemes. Such methods have excellent conservation properties, can be readily designed for arbitrary order of accuracy, and have the potential for implementation on unstructured meshes with adaptivity. To ensure the positivity of the solution, one can use a maximumprinciplesatisfying limiter that has been recently proposed by Zhang and Shu in [52] for conservation laws on cartesian meshes, and later extended on triangular meshes [56]. This limiter has been used to develop positivitypreserving schemes for compressible Euler [53, 55], shallow water equations [48], and VlasovBoltzmann transport equations [10]. It has also been employed recently in the framework of semiLagrangian DG methods [43, 41] for the VP system.
The scope of the present paper is as follows: we focus on a detailed study of the RKDG scheme for the Vlasov equation from both the numerical and analytical points of view. Since we are only considering onedimensional problems, we use the classical representation of the solution by Green’s function to compute the Poisson equation; therefore, the electric field is explicitly given as a function of the numerical density. This removes all discretization error from the Poisson equation and lets us more accurately investigate our DG solver for the Vlasov equations. We rigorously study recurrence, which is an important numerical phenomenon that commonly appears with many solvers. We use Fourier analysis and obtain eigenvalues of the amplification matrix, and then investigate the impact of different polynomial spaces on the quality of the solution by examining conserved quantities as well as convergence to BGK states [5] for some choices of initial states. We consider benchmark test problems such as simulations of Landau damping phenomena for the linearized and nonlinear Vlasov Poisson systems, twostream instability, and their long time BGK states and the formation of KEEN waves, both for the nonlinear system as well.
The remaining part of the paper is organized as follows: in Section 2, we describe the numerical algorithm and summarize its conservation properties. In Section 3, we study the recurrence phenomena that occurs for linear Vlasov type transport equations discretized by DG methods with various polynomial orders. Sections 4.1 and 4.2 are devoted to discussions of simulation results for the linearized and nonlinear VP system, respectively, for diverse choices of initial data and external drive potentials. We conclude with a few remarks in Section 5.
2 Numerical methods
In this section we first describe the proposed DG numerical algorithm and then discuss some of its basic conservation properties related to the quantities of (1). This is done for both the fully nonlinear VP system of (1) and the linearized VP system obtained by linearizing about the homogeneous equilibrium , with corresponding vanishing electric equilibrium field. Periodic boundary conditions in are assumed.
Thus, setting and expanding the system to first order approximation, the perturbation satisfies the Linear VlasovPoisson (LVP) system,
(3)  
where has been replaced by to ease the notation. We find it convenient and efficient to intertwine the discussion of our algorithms for the VP and LVP systems. To avoid confusion in Section 2.1 we underline the words linear and nonlinear, signaling where discussions specific to each apply.
2.1 Numerical algorithm
For onedimensional problems, we use a mesh that is a tensor product of grids in the and directions, because this simplifies the definitions of the mesh and polynomial space for the Poisson equation. Specifically, the domain is partitioned as follows:
where is chosen appropriately large to guarantee for . This is a reasonable assumption, because of the wellposedness of the onedimensional VlasovPoisson system as indicated in [21]. The grid is defined as
where and are center points of the cells.
We will make use of several approximation spaces. For the domain, we consider the piecewise polynomial space of functions ,
where is the space of polynomials in one dimension of degree up to . For the space, we consider two approximation spaces of functions ,
and
where denotes all polynomials of degree at most in and on , and . These two spaces are widely considered in the DG literature for multidimensional problems. A simple calculation shows that the number of degrees of freedom of is . For this is larger than the number of degrees of freedom of , which is .
First we describe the RKDG scheme for the linear Vlasov equation. We seek (or ), such that
(4) 
holds for any test function (or ). Here and below, we use the following notations: is the discrete electric field, which is to be computed from Poisson’s equation, , , and is a numerical flux. We can assume that in each , has a single sign by properly partitioning the mesh. Then, the upwind flux is defined as
(5) 
The scheme for the nonlinear Vlasov equation is similar. We seek (or ), such that
(6) 
holds for any test function (or ). The upwind flux for has been defined in (5) and the new flux needed for the nonlinear case is given by
(7) 
The above descriptions coupled with a suitable time discretization scheme, such as the TVD RungeKutta method [44], completes the RKDG methods. For example, suppose the semidiscrete schemes in (2.1) and (2.1) are written in the compact form
where for the linear Vlasov of (2.1)
while for the nonlinear Vlasov of (2.1)
The third order TVDRK method implements the following procedure for going from to :
(8)  
Poisson’s equation is used to obtain , , and . Beyond periodicity, we need to enforce some additional conditions to uniquely determine . For example, we can set one end of the spatial domain to ground, i.e. set . In the onedimensional case, then the exact solution can be obtained if we enforce . For the nonlinear system we obtain
where , and
(9) 
while for the linear system Poisson’s equation gives
where , and
(10) 
From (9) and (10), we see that if (or ), then ; hence, . The above approach uses the classical representation of the solution by Green’s function and will be referred to as the “exact” Poisson solver. It is valid only for the onedimensional case. For higher dimensions, a suitable elliptic solver needs to be implemented, such as those discussed in [23]. Here we use the exact solver to entirely eliminate discretization error from Poisson’s equation and, thereby, spotlight the performance of the Vlasov solver.
Below we describe positivitypreserving limiters, as summarized in [54]. We only use such a limiter to enforce the positivity of for the nonlinear VP system. For each of the forward Euler steps of the RungeKutta time discretization, the following procedures are performed:

On each cell , we evaluate , where , and denote the Gauss quadrature points, while denote the GaussLobatto quadrature points.

We compute , where is the cell average of on , and . This limiter has the effect of maintaining the cell average, while “squeezing” the function to be positive at all points in .

Finally, we use instead of to compute the Euler forward step.
This completes the description of the numerical algorithm.
2.2 Scheme Conservation properties
In the following, we will briefly review and discuss some of the conservation properties of the above RKDG scheme for the nonlinear VP equations without the positivitypreserving limiter. Some of those results have been reported in [22] and [3].
Proposition 1
(charge conservation) For both the and spaces,
which implies
for the fully discrete scheme (2.1). Here,
denotes the contribution from the phase space boundaries located at , and should be negligible if is chosen large enough.
Remark: Charge conservation (or mass conservation or probability normalization as it is sometimes called) states that the total charge will be preserved on the discrete level up to approximation errors associated with the phase space boundaries. The proof is straightforward and, therefore, omitted. The same conclusion can be proven for the linearized system. The positivity preserving limiter does not destroy this property because it keeps the cell averages unchanged.
Proposition 2
(Semidiscrete stability – enstrophy decay) For both the and spaces, . Hence,
The proof, for an arbitrary field , can be found in [10], Theorem 4, which applies directly here by setting the collisional form in that proof.
For the remainder of this subsection we will assume the DG solution satisfies the velocity boundary conditions . This is a reasonable assumption when is large enough. In particular, this will guarantee exact charge conservation, which implies that is constant in time . Therefore, using the definition of in (9), we can obtain periodicity in , i.e, . Without this assumption the propositions below contain multiple boundary terms and the proof becomes technical.
Proposition 3
(Momentum conservation) Assuming , for both the and spaces when , , which implies
for the fully discrete scheme.
Proof. Choosing in (2.1), we have
and using the exact Poisson solver together with the periodicity of and yields the following:
which completes the proof.
Remark: The above proof holds for the linearized system as well. Note, however, it relies on the use of the exact Poisson solver. For a full numerical DG Poisson solver, such as that developed in [23] for the discretization Poisson equation, exact momentum conservation remains true, as was proven in [22] by means of the DFUG method developed there for dealing with the discretized Poisson equation. However, the positivitypreserving limiter we use here will destroy this property because it was not designed to conserve the numerical momentum.
Proposition 4
(Semidiscrete total energy equality) Assuming , for both the and spaces when ,
where the operator and is any projection such that and at , for .
Proof. Choosing in (2.1) yields
and
where in the second line, we have used the periodicity and continuity of and . Therefore, we have proven that
On the other hand, upon choosing in (2.1) and using the periodicity and continuity of , we can verify that . The exact solution obviously satisfies from the continuity and periodicity of , and therefore we are done.
The above proof indicates that the variation in the total energy will be related to the error of the solution, , together with the projection error, . In [22, 23], error estimates for DG schemes with NIPG methods for the Poisson equations are provided for multiple dimensions. In [3], optimal accuracy of order for the semidiscrete scheme with spaces has been proven under certain regularity assumptions. We remark that in [3] conservation of the total numerical energy is proven when the Poisson equation is solved by a local DG method with a special flux. Unfortunately, no numerical simulations of linear Landau damping or of the nonlinear VP system, such as those done in [23] or in Section 4 of this present manuscript, have been performed up to this date by the scheme proposed in [3], so a comparison is not possible.
3 On recurrence
In this section we study recurrence, a numerical phenomenon that is known to occur in simulations of Vlasovlike equations. Its study is important because it provides information about the numerical accuracy of a scheme. Recurrence was observed in the semiLagrangian scheme of Cheng and Knorr [8], where a simple argument for its occurrence was provided. In this section, we carry out a detailed study of recurrence for the DG method.
We study recurrence of our algorithm applied to the linear advection equation on , since it is tractable and contains the basic recurrence mechanism; results for the full Vlasov system tend to be quite similar. The initial condition we consider is , and the equilibrium distribution is taken to be either the Maxwellian or Lorentzian distribution, viz.
For the Maxwellian equilibrium, , we take , and for the Lorentzian equilibrium, , we take .
The exact solution for the advection equation is . Hence, a simple calculation shows for the Maxwellian distribution; similarly, for the Lorentzian, . Thus, we see how the density for each case should decay to zero. The failure of decay and the occurrence of recurrence noted for the semiLagrangian scheme of [8] stems from the finite resolution in the velocity space and, indeed, the recurrence time depends on the mesh size in .
To be specific, we repeat the definition of DG scheme for this equation, which amounts to (2.1) with set to zero: we find (or ) , such that
(11) 
holds for any test function (or ). Again is the upwind numerical flux of (5). In the analysis below, we always assume time to be continuous, because recurrence is mainly a phenomenon that comes from the spatial and velocity discretization.
3.1 The case of
For the piecewise constant case, the DG method is equivalent to a simple first order finite volume scheme and we can derive rigorously the behavior for . Suppose we define on cell , and assume uniform grids in both directions. Moreover, we assume to be even for simplicity. With this assumption, the location of the cell center is . Now (11) simply becomes
(12) 
The initial condition chosen is clearly equivalent to . We prove that the scheme above gives
(13) 
where is given in (14) below.
Upon plugging (13) into (12), we have
Hence,
which can be summarized as
(14) 
Therefore, the real part of is always negative, this means the magnitude of will always damp, but because of the dependence it does so at different rates for different cells. Since the density
the density will damp at a rate between and . Another important fact is that recurring local maxima of the density will have a period that satisfies . If we define , then . When , , and this coincides with the recurrence time obtained in [8].
Next we compare the above theoretical prediction with numerical results. In all of the calculations below, we take , and the mesh size to be . In Figure 1, we display results for numerical runs using piecewise constant polynomials and time discretization using TVDRK3. (We use the third order method to minimize the time discretization error.) We plot in Figure 1. First, we notice the pattern of has the expected periodic structure with damping for both Maxwellian and Lorentzian equilibria. For the Maxwellian distribution, a simple calculation yields . Similarly, with the formulas above, the smallest damping rate is , while the biggest is . For the Lorentzian distribution,