# A discontinuous Galerkin method for the Vlasov-Poisson system

## Abstract

A discontinuous Galerkin method for approximating the Vlasov-Poisson system of equations describing the time evolution of a collisionless plasma is proposed. The method is mass conservative and, in the case that piecewise constant functions are used as a basis, the method preserves the positivity of the electron distribution function and weakly enforces continuity of the electric field through mesh interfaces and boundary conditions. The performance of the method is investigated by computing several examples and error estimates associated system’s approximation are stated. In particular, computed results are benchmarked against established theoretical results for linear advection and the phenomenon of linear Landau damping for both the Maxwell and Lorentz distributions. Moreover, two nonlinear problems are considered: nonlinear Landau damping and a version of the two-stream instability are computed. For the latter, fine scale details of the resulting long-time BGK-like state are presented. Conservation laws are examined and various comparisons to theory are made. The results obtained demonstrate that the discontinuous Galerkin method is a viable option for integrating the Vlasov-Poisson system.

## 1Introduction

The single species Vlasov-Poisson system is a nonlinear kinetic system that models time evolution of a collisionless plasma consisting of electrons and a uniform background of fixed ions under the effects of a self-consistent electrostatic field and possibly an externally supplied field. The Vlasov equation models the transport of the electrons that are coupled to the electrostatic potential through Poisson’s equation.

The unknown electron distribution function, a phase space density, is denoted by where the independent variables are position, velocity, and time, respectively. For a given , the quantity denotes the number of electrons contained in the infinitesimal phase space volume element centered about at time . Upon a proper renormalization, can be interpreted as a probability distribution function for the electrons over the phase space.

The Vlasov-Poisson system has solutions that can exhibit a variety of dynamical phenomena [69]. One of the most well-known effects is filamentation or ‘phase mixing’ as it is sometimes called, which occurs when different characteristics surfaces associated to the nonlinear transport (Vlasov) equation wrap in the phase space. This effect results in stiff gradients, since generally takes on disparate values along different characteristics.

Related to filamentation is the interesting phenomenon of Landau damping [42]: electrostatic disturbances can be interpreted as the interaction between plasma waves and electrons with a resulting net energy transfer from the wave to electrons leading to an exponential collisionless damping of the electric field modes in time. Because of such phenomena, the Vlasov-Poisson equation can be challenging to simulate numerically.

Existing numerical techniques for solving the Vlasov-Poisson system can be divided into two groups: (i) those that approximate the system in the phase space directly and (ii) those that transform the system into a different coordinate space. The numerical approaches that treat the phase space directly do not, however, usually involve discretizing the Vlasov equation. Rather most of these methods take advantage of the characteristic structure of the Vlasov equation, which implies that the plasma particles evolve along trajectories that satisfy a given set of ordinary differential equations. The most widely used particle scheme is the Particle-in-Cell (PIC) method [9], which represents ensembles micro-particles as a finite number of macro-particles. Each macro-particle is then assumed to evolve along a characteristic trajectory, where the electric field defining the trajectory is computed via any standard scheme. The PIC method seems to give reasonable results in cases where the tail of the distribution is negligible and a large number of particles is not necessary. Otherwise, the method suffers from numerical noise that is proportional to where is the number of particles.

Other methods based on the discretization of the phase space have also been proposed. In [11], an operator splitting method was introduced and shown to be both efficient and accurate for solving a wide range of problems. A continuous finite element method was developed in [72] and was shown to achieve results similar to those obtained in [11]. A positive and mass conservative scheme was employed in [33] to solve both linear and nonlinear damping problems in one- and two-dimensional physical space. This method is defined at a given time step by first building a piecewise constant approximation over a mesh of the phase space using the approximation obtained from the previous time step and two correction terms whose values are found by solving two fixed point problems. The piecewise approximation is then used in conjunction with a slope limiter to reconstruct a local polynomial approximation of for each cell in the mesh.

Transform methods based on Fourier or Hermite series have also been used, e.g. in [36] and more recently in [31]. In [40], the phase space was transformed using a Fourier transform and a splitting method was employed to advance the approximation in time. This method also included a filamentation filtering step for the purpose of smoothening the filaments. The numerical results obtained using this method seem reasonable only for problems in which filamentation is not a dominant effect, where perhaps the nonlinearity either slows down or prevents the onset of filamentation. However, this method may be inadequate for problems where the physics of interest depends upon the filamentary nature of the distribution, such as is the case for Landau damping.

The objective of this paper is to propose a coupled Upwind Penalty Galerkin (UPG) method for the approximation of the Vlasov-Poisson system and to evaluate its numerical efficacy. Our UPG method gives a unified approach for approximating both the hyperbolic and elliptic parts of the Vlasov-Poisson system. Specifically, the Vlasov equation is discretized using the standard Upwind Galerkin (UG) scheme for conservation laws [24] and the Poisson equation is discretized using one of the three Discontinuous Galerkin (DG) interior penalty available schemes [70]. Stability and convergence estimates for the UPG method were presented in [37] for the six-dimensional phase space In the same reference, the method was shown to be both locally and globally mass conservative.

More specifically, the semi-discrete UPG approximation to the electron distribution function is defined to be the solution of a first-order, nonlinear, ordinary differential equation (ODE) system. Moreover, it has been shown that the method preserves positivity of when piecewise constants basis functions are used to approximate the solution to the Vlasov-Poisson system.

In this manuscript we show the numerical efficacy of the DG method by performing accuracy, convergence, and conservation tests on computed UPG approximate solutions to a variety of linear and nonlinear problems where sufficiently data or information of the true nonlinear solutions has been established. The computed results for these problems are benchmarked against known theoretical results and are compared to results obtained using established methods.

For computing plasma problems the UPG method offers several advantages. In particular, the local nature of the method facilitates adaptive mesh refinements with easy adaptation to parallelization techniques. By taking advantage of these benefits, regions of the phase space where the electron distribution experiences strong filamentation or boundary layer effects can be resolved by local mesh refinements. The discontinuous nature of the method also helps to resolve the stiff gradients associated with filamentation, since requiring the approximation to be continuous in these cases can be too restrictive and typically lead to excessive numerical diffusion and oscillatory behavior. Due to the fact that the method imposes boundary conditions weakly, a variety of boundary conditions can easily be accommodated.

Recently, an alternative DG formulation for the Boltzmann-Poisson system was introduced in [12] for simulations of hot electron transport for one and two space dimensional nanometer scale devices, where kinetic corrections are known to be very significant. In [20] a DG scheme was constructed for the Vlasov-Boltzmann equation by means of a maximum principle satisfying a limiter for conservation laws [74], and it was shown to be high order accurate and positivity-preserving, not only for piecewise constant basis functions but also for higher order polynomial approximations. Consequently a new approximation of the Vlasov system based on a UPG scheme for higher order polynomial basis functions is currently being developed and tested [18]. Thus, the DG method is well-suited to approximate a range of plasmas spanning from the collisionless to the highly collisional regimes.

For completeness of this introduction we include some historical remarks regarding DG methods. The initial development of these numerical approximating schemes for hyperbolic and elliptic equations occurred independently, but nearly simultaneously. In 1973, the first DG scheme for linear hyperbolic equations was introduced by Reed and Hill for approximating a neutron transport equation [57]. This work was followed by Lesaint and Raviart [43] in 1974, where *a priori* error estimates were proved for the DG method applied to two-dimensional, linear hyperbolic problems. The DG schemes for hyperbolic problems were further studied in the series of papers [24], which culminated in the introduction of the local discontinuous Galerkin (LDG) method [25]. The generality of the LDG method was further extended to the multidimensional setting under more relaxed assumptions on the data in [21]. One of the first DG schemes for approximating the solutions to second-order elliptic equations was introduced in 1971 by Nitsche [55] where Dirichlet boundary conditions were enforced weakly rather than strongly through the use of a penalty term. Shortly thereafter, applications of the penalty method to Laplace’s equation were proposed by Babuška et al. in [4].

The use of penalty terms across interior faces as a means of enforcing continuity among adjacent elements was introduced in [70] and [71] using a symmetric interior penalty Galerkin (SIPG) finite element method. A non-symmetric interior penalty method (NIPG) similar to the SIPG method was proposed and analyzed in [59]. The incomplete interior penalty method (IIPG) was introduced in [63] and is very similar to the SIPG and NIPG methods.

The outline of this paper is as follows. In Section 2, we describe the Vlasov-Poisson system and give a brief discussion of Landau damping. In Section 3, the UPG method for the approximation of the Vlasov-Poisson system is derived and the error estimates associated to the approximation of the system are stated. In Section 4, several numerical simulation results are presented and analyzed, including a study of the free streaming operator (simple advection), Landau damping for Maxwellian and Lorenztian equilibria, strong nonlinear Landau damping and a careful study of a symmetric two-stream instability, all for the case of repulsive potential forces. In Section 5, we comment on the efficiency of the UPG method, draw some conclusions, and remark on future work. Finally, an Appendix dedicated to the analysis of dispersion relations for the Lorentzian and two-stream equilibria is included.

## 2The Vlasov-Poisson System

The Vlasov-Poisson system considered in this work has been scaled by the usual characteristic time and length scales, i.e., time is scaled by the inverse plasma frequency , length by the Debye length , and velocity accordingly by a thermal velocity .

Using this nondimensionalization, the Vlasov-Poisson problem is as follows: For the divergence free vector in phase space

find the electron distribution function and the electric field with corresponding electrostatic potential such that, for fixed ,

subject to an initial condition on and boundary conditions on and The domain of definition of the initial boundary value problem , where the physical domain can either be bounded or all of with . The boundary condition given for depends on If then the condition must be enforced both as and as If is bounded, then a condition must be imposed on along the inflow boundary defined by

with being the unit outward normal vector to Often is given in some parts of the boundary and may be periodic in other parts of the boundary region . In this manuscript we assume periodic boundary conditions in space and the decaying boundary condition in velocity.

The Poisson equation must also be endowed with spatial boundary conditions, either Dirichlet, Neumann, Robin, or periodic, on different disjoint regions of the boundary . We denote the Dirichlet portion of the boundary by . If the measure of the Dirichlet boundary is zero, i.e. , and one has homogeneous or periodic boundary conditions such that , then in order to maintain a well-posed problem that keeps the existence and uniqueness of the corresponding Poisson boundary value problem, one needs to add to the solution space the compatibility (neutrality) condition , or equivalently on each connected component of the spatial domain .

Macroscopic fluid quantities of interest are easily computed from . The electron density , current density , kinetic energy density and electrostatic energy density are defined by

These quantities satisfy a number of respective conservation laws (see e.g. [58]). In particular, it is well-known that the Vlasov-Poisson system conserves total particle number, momentum, energy, and the Casimir invariants, which are given, respectively, by

The notation in ( ?) refers to an arbitrary function of and includes the ‘enstrophy’ when , entropy when , or particle number, as in (Equation 4), when . We will check the invariance of these quantities in our nonlinear computations, particularly in Section ?.

Interesting properties of the Vlasov-Poisson system result by considering a linear perturbation to an equilibrium distribution over the 2D-phase space Specifically, suppose that , where is a given equilibrium probability distribution, and are -periodic in , and the initial average value of over is zero. Equations ( ?)-( ?) imply that satisfies

Supposing and dropping this term from ( ?) leads to

The linear system of ( ?)-( ?) was analyzed by using the Laplace transform in the famous paper of Landau [42], by expansion in terms of continuum eigenfunctions in [69], and by a tailored integral transform introduced in [51]. Landau showed that an electric field mode decays exponentially in the long-time limit, which we investigate in Section ? for two well-known equilibria: the Maxwellian

and the Lorentzian,

## 3Method Formulation

In this section, we derive the UPG method for the Vlasov-Poisson system. The derivation proceeds by first discretizing the Vlasov equation using the standard UG discretization for transport equations [21]. Here it is assumed that the electric field is given and hence the divergence-free flow field , defined in for the Vlasov equation , is known. Afterwards, the DG discretization for the Poisson equation is considered.

It is assumed that any mesh for the phase space, where by mesh we mean a partitioning of the phase space into convex sets called elements, is the Cartesian cross product of a mesh for the physical domain and a mesh for the velocity domain. Under this assumption, the physical and velocity domains can be independently refined. Given a mesh for the phase space, the UG method then defines an approximate solution to the true Vlasov solution in such a way that at any given time the approximate solution restricted to each element of the mesh is a polynomial function. However, the approximate solution is not required to be continuous across the intersections of any two adjacent elements, so it is a piecewise defined polynomial function with respect to the mesh at any given time.

In order to compute the approximation to the electrostatic potential from Poisson’s equation , for a given distribution function, one may use one of three interior penalty methods that weakly enforce both approximate continuity across the interior mesh faces and Dirichlet boundary regions. These three alternative methods, symmetric interior penalty Galerkin (SIPG) [70], non-symmetric interior penalty Galerkin (NIPG) [59], and incomplete interior penalty Galerkin (IIPG) [63] are discussed in detail and *a priori* error estimates for each of them are given in the respective references. The only difference among the three methods is in the value of one specific parameter that arises in the weak formulation that is common to each of them. Thus, for a given space charge function, each penalty Galerkin method defines a piecewise polynomial approximation of the true solution to the Poisson equation using a mesh in the spatial domain.

The spatial domain mesh used in the discretization of Poisson’s equation is required to be the same as that used in the UG discretization of the Vlasov equation. This requirement is a practical one, both in terms of analysis and implementation. However, the polynomial degree of the potential approximation on a given element of the spatial mesh is not required to equal the degree, with respect to , of the polynomial approximation of the distribution .

Thus, the UPG method of approximation to the Vlasov-Poisson system is defined by coupling together the UG method of approximation to the Vlasov equation with interior penalty methods of approximation to the Poisson equation. This nonlinear semi-discrete approximation results in a first-order nonlinear ODE system, the solution of which determines the approximation of . The resulting ODE system is readily solved using an explicit conservative time-integrator such as the Runge-Kutta method. Moreover, in the process of computing , both the approximation of the electric field and the approximation of the potential are computed by one of the penalty methods for the elliptic equation.

### 3.1Preliminaries

We assume the computational domain in velocity space is a bounded set and that the approximate solution to is assumed to vanish in . It is then, implicitly for our simulation, assumed that the velocity support of the approximation to the true solution of the Vlasov-Poisson system is contained in for all times. This is a reasonable assumption for problems with spatial periodic boundary conditions, as it is expected that most of the density associated with the approximation to will be contained in a sufficiently large fixed set . The error due to this assumption only depends on the density of the true solution computed on the complement of in .

Further, the conservative nature of the transport equation with the space dependent divergence-free flow field , motivates us to choose the computational scheme as follows:

Let be a sequence of successively refined meshes of the bounded domain and let be a sequence of successively refined meshes of the bounded domain , where Given the meshes and , the elements and comprising each of the respective meshes are sets of the following types: intervals, if ; triangles or quadrilaterals, if ; and tetrahedra, prisms or hexahedra, if . The corresponding spatial refinement level and velocity refinement level are defined by and , respectively.

A sequence of successively refined meshes of the, now, computational domain is generated by defining each mesh to be where the refinement level is . Thus, for any given element there exists a unique pair of elements and such that , which is equivalent to the existence of an invertible mapping from to , where .

The derivation of the following UPG method requires the use of the broken Sobolev space , , which is defined as follows:

i.e., is the space of those functions that have elementwise weak derivatives up to, and including, the order . Then for nonnegative integers and , the discontinuous approximation space is given by

where denotes the space of polynomials on a set with degree less than or equal to in each variable. Thus, , where denotes the space of polynomials satisfying that the sum of the degrees of all the variables is less than or equal to .

The choice of for basis functions is suitable for Cartesian meshes in both -space and -space, respectively, where trace and inverse inequalities that are derived by mapping to the reference element in the approximating framework are possible, as first introduced in [37].

However, one may use triangles in two-dimensions, and prisms or hexahedra in three-dimensions, for both -space and -space, for which the natural choice of polynomial space would be . We point out that these selections of approximating spaces is consistent with the divergence free, linear, and conservation form of the Vlasov equation, and the fact that the choice of the mesh associated with the computational domain is a set of product mesh elements , for , and . Such is clearly preserved by the Vlasov flow and the corresponding approximation results for will be valid. This choice may be preferable in higher dimensions since the number of degrees of freedom of the basis functions in is , while for is , for -dimensional calculations.

The discontinuous nature of the space needs the introduction of mesh faces. If is a boundary element, then is called a boundary mesh face. If and are two intersecting elements whose common intersection lies in the interior of , then is said to be an interior mesh face. The set of all mesh faces is denoted by , where is an interior face if and a boundary face if . Each face is associated with a unit normal vector . For , is chosen to be the outward unit normal to . For , we fix to be one of the two unit normal vectors to For every interior face the elements and will always be used to denote the two unique elements such that Moreover, it is always assumed that satisfies on , where denotes the outward unit normal to Then on .

The fact that each mesh is the product of a spatial mesh and a velocity mesh gives a specific structure to the boundaries of the elements and to the set of mesh faces . It follows that for we have . We denote the set of mesh faces for and by and , respectively, where is an interior face if and a boundary face if , and is an interior face if and a boundary face if . Then, given any arbitrary , either there exist an and a such that or there exist a and an such that .

When considering functions in , we use the usual average and jump operators, respectively, defined for along an interior face by

The above definitions are also valid for vector-valued functions , in which case it follows that .

### 3.2Upwind Galerkin Approximation of the Vlasov Equation

Here we describe the UG scheme for the Vlasov equation in full generality for a -dimensional ( or ) phase space with inflow boundary conditions and piecewise polynomials of arbitrary degree approximating . The simpler derivation of the method for a two-dimensional phase space with periodic boundary conditions in and a piecewise constant approximation to is given explicitly in Section 3.5, since this method was used for the numerics of Section 4. Note, for this simpler version the derivation does not require use of the set of mesh faces .

For a given time and data trio , the Vlasov equation along with the corresponding initial and boundary conditions are

where , is assumed given, , and . It is assumed that and where are fixed. Following , we also define (keeping the same notation without loss of generality) the computational inflow boundary , associated with the computational domain , by

where is the outward unit normal to . Then, .

For domains , let denote the -inner product. To distinguish integration over domains , we use the notation . A weak formulation for ( ?)-( ?) is derived by multiplying equation ( ?) by an arbitrary test function and integrating by parts over an arbitrary . This yields

with being the outward unit normal to and denoting the interior trace of , i.e., .

Upon summing equation (Equation 10) over all and weakly enforcing the inflow boundary condition, we get

Approximating by a function , which may be discontinuous across the interior faces, requires the introduction of a numerical flux . A standard technique is to replace by its “upwind value” on the interior faces [21], where along each interior face is defined by the given advecting flow field according to

Consistency follows from condition (Equation 12), since on whenever is continuous across each interior face . We also note that depends nonlinearly on as, in general, if and are two flow fields having different values on and if is discontinuous across , then . Replacing on the interior faces in (Equation 11) by its upwind value leads to

which is the UG scheme for the Vlasov equation. Finally, defining the bilinear operator by

and the corresponding linear operator , depending on the inflow data , by

yields a variational formulation for the semi-discrete problem of finding the approximation to , satisfying,

for all , with and approximations of the initial data and the inflow boundary data on , respectively.

We note that (Equation 16) produces a first-order ODE system to be described below. Indeed, for each let be a basis for , where and then extend the domain of these functions to by defining each to be identically zero in . If denotes the unique vector such that the UG approximation satisfies

then . Inserting (Equation 17) into (Equation 16) yields

Finally, since is a basis for , then (Equation 18) is equivalent to

where contains the indices of all neighboring elements of .

Equation (Equation 19) is seen to generate an equivalent matrix system, where rows of the matrix are generated at a time by sequentially taking to equal and for each sequentially taking to equal in Eq. (Equation 19). This procedure results in the matrix ODE system

where is a constant matrix and is the corresponding sparse matrix, both of which are of dimension , and is a vector of length .

Since the support of the functions is , , it follows that is a block-diagonal matrix, where each block is an matrix. This means that the inverse of is easily computed. Thus, the UG approximation is equivalently defined to be the unique solution to

where the initial condition is uniquely determined by ( ?). To solve this system in time, a conservative explicit time stepping method such as the Runge-Kutta method can be used.

### 3.3Interior Penalty Approximations of the Poisson Equation

In order to make this manuscript self contained, we also discuss the interior penalty formulations for the Poisson equation that weakly enforces approximate continuity across interior mesh faces and Dirichlet boundary regions. As already noted, the mesh used to discretize the Poisson equation must be the same as that for the spatial domain used for the Vlasov equation, but the polynomial degree for the Poisson equation need not be equal to the degree in for .

Hence, for a given i) source , ii) boundary data in the portion of the boundary referred as the Dirichlet boundary , and iii) (homogeneous Neumann) or periodic boundary conditions on , the boundary complement of the Dirichlet region, then the more general form of Poisson equation with a positive definite permittivity term discretized using the symmetric interior penalty (SIPG) method [70], incomplete interior penalty (IIPG) method [63], or non-symmetric interior penalty (NIPG) [59] method is

where is any positive-definite continuous function in .

For , =1,2, or 3, recall is the -inner product with integration over , while the notation is for boundary integrals. In addition, we use the identity

where is the outward unit normal to .

Thus, the corresponding schemes SIPG, IIPG, and NIPG are all derived by multiplying ( ?) by an arbitrary test function , , integrating by parts on each and summing each of the resulting local equations. Whence, one obtains the non-symmetric variational formulation given by

where the identity (Equation 22) was used to represent the inner boundary integrals.

In particular, since the true solution for a bounded right-hand-side has at least the regularity , it follows that and along every interior face [32]. Thus, In order to get a good approximation for a regular solution satisfying these two jump conditions, one adds ‘interior’ penalty terms that vanishes on the true solution . These terms are of the form

where the values selected for for the SIPG, IIPG and NIPG methods are -1, 0 and 1, respectively. Similarly, such penalization is also required for the Dirichlet boundary region as follows.

Indeed, we add to the bilinear form the interior penalty terms

which are now completed by weakly enforcing both approximate continuity across the interior mesh faces and Dirichlet boundary regions , by including the penalization

where is an arbitrary penalty parameter that is usually set equal to unity. Note that homogeneous or periodic boundary conditions vanish on boundary terms of the corresponding bilinear structure. Consequently, they do not require the boundary penalization term.

Usually the penalization terms are denoted by the following non-symmetric bilinear form

Adding both penalty terms and (Equation 24) to the left-hand-side of (Equation 23) results in

Equation (Equation 28) completely defines each of the three interior penalty schemes.

Setting equal to the first four terms on the left-hand-side of (Equation 28), where dependence on the parameter from is noted as a subscript, we let the bilinear operator and the linear operator be equal to the two terms on the right-hand-side of (Equation 28). Then the function is the corresponding interior penalty Galerkin approximation to the Poisson solution , if