Discontinuous Galerkin Methods for the Vlasov-Maxwell Equations

# Discontinuous Galerkin Methods for the Vlasov-Maxwell Equations

## Abstract

Discontinuous Galerkin methods are developed for solving the Vlasov-Maxwell system, methods that are designed to be systematically as accurate as one wants with provable conservation of mass and possibly total energy. Such properties in general are hard to achieve within other numerical method frameworks for simulating the Vlasov-Maxwell system. The proposed scheme employs discontinuous Galerkin discretizations for both the Vlasov and the Maxwell equations, resulting in a consistent description of the distribution function and electromagnetic fields. It is proven, up to some boundary effects, that charge is conserved and the total energy can be preserved with suitable choices of the numerical flux for the Maxwell equations and the underlying approximation spaces. Error estimates are established for several flux choices. The scheme is tested on the streaming Weibel instability: the order of accuracy and conservation properties of the proposed method are verified.

V

lasov-Maxwell system, discontinuous Galerkin methods, energy conservation, error estimates, Weibel instability

{AMS}

65M60, 74S05

## 1 Introduction

In this paper, we consider the Vlasov-Maxwell (VM) system, the most important equation for the modeling of collisionless magnetized plasmas. In particular, we study the evolution of a single species of nonrelativistic electrons under the self-consistent electromagnetic field while the ions are treated as uniform fixed background. Under the scaling of the characteristic time by the inverse of the plasma frequency , length by the Debye length , and electric and magnetic fields by (with the electron mass, the speed of light, and the electron charge), the dimensionless form of the VM system is

 ∂tf +ξ⋅∇xf+(E+ξ×B)⋅∇ξf=0 , (1.1a) ∂E∂t =∇x×B−J,∂B∂t=−∇x×E , (1.1b) ∇x⋅E =ρ−ρi,∇x⋅B=0 , (1.1c)

with

 ρ(x,t)=∫Ωξf(x,ξ,t)dξ,J(x,t)=∫Ωξf(x,ξ,t)ξdξ ,

where the equations are defined on , denotes position in physical space, and in velocity space. Here is the distribution function of electrons at position with velocity at time , is the electric field, is the magnetic field, is the electron charge density, and is the current density. The charge density of background ions is denoted by , which is chosen to satisfy total charge neutrality, . Periodic boundary conditions in -space are assumed and the initial conditions are denoted by , and . We also assume that the initial distribution function , i.e., is in a Sobolev space of order and is integrable with finite energy in -space, where . The initial fields and are also assumed to be in .

The VM system has wide importance in plasma physics for describing space and laboratory plasmas, with application to fusion devices, high-power microwave generators, and large scale particle accelerators. The computation of the initial boundary value problem associated to the VM system is quite challenging, due to the high-dimensionality (6D+time) of the Vlasov equation, multiple temporal and spatial scales associated with various physical phenomena, nonlinearity, and the conservation of physical quantities due to the Hamiltonian structure [43, 44] of the system. Particle-in-cell (PIC) methods [5, 35] have long been very popular numerical tools, in which the particles are advanced in a Lagrangian framework, while the field equations are solved on a mesh. This remains an active area of research [22]. In recent years, there has been growing interest in computing the Vlasov equation in a deterministi c framework. In the context of the Vlasov-Poisson system, semi-Lagrangian methods [11, 55], finite volume (flux balance) methods [6, 23, 24], Fourier-Fourier spectral methods [39, 40], and continuous finite element methods [58, 59] have been proposed, among many others. In the context of VM simulations, Califano et al. have used a semi-Lagrangian approach to compute the streaming Weibel (SW) instability [9], current filamentation instability [42], magnetic vortices [8], magnetic reconnection [7]. Also, various methods have been proposed for the relativistic VM system [54, 4, 56, 36].

In this paper, we propose the use of discontinuous Galerkin (DG) methods for solving the VM system. What motivates us to choose DG methods, besides their many widely recognized desirable properties, is that they can be designed systematically to be as accurate as one wants, meanwhile with provable conservation of mass and possibly also the total energy. This is in general hard to achieve within other numerical method frameworks for simulating the VM system. The proposed scheme employs DG discretizations for both the Vlasov and the Maxwell equations, resulting in a consistent description of the distribution function and electromagnetic fields. We will show that up to some boundary effects, depending on the size of the computational domain, the total charge (mass) is conserved and the total energy can be preserved with a suitable choice of the numerical flux for the Maxwell equations and underlying approximation spaces. Error estimates are further established for several flux choices. The DG scheme can be implemented on both structured and unstructured meshes with provable accuracy and stability for many linear and nonlinear problems, it is advantageous in long time wave-like simulations because it has low dispersive and dissipative errors [1], and it is very suitable for adaptive and parallel implementations. The original DG method was introduced by Reed and Hill [51] for a neutron transport equation. Lesaint and Raviart [41] performed the first error estimates for the original DG method, while Cockburn and Shu in a series of papers [18, 17, 16, 15, 19] developed the Runge-Kutta DG (RKDG) methods for hyperbolic equations. RKDG methods have been used to simulate the Vlasov-Poisson system in plasmas [34, 33, 13] and for a gravitational infinite homogeneous stellar system [12]. Some theoretical aspects about stability, accuracy and conservation of these methods in their semi-discrete form are discussed in [33, 3, 2]. Recently, semi-Lagrangian DG methods [52, 50] were proposed for the Vlasov-Poisson system. In [37, 38], DG discretizations for Maxwell’s equations were coupled with PIC methods to solve the VM system. In a recent work [?], error estimates of fully discrete RKDG methods are studied for the VM system.

The rest of the paper is organized as follows: in Section 2, we describe the numerical algorithm. In Section 3, conservation and the stability are established for the method. In Section 4, we provide the error estimates of the scheme. Section 5 is devoted to discussion of simulation results. We conclude with a few remarks in Section 6.

## 2 Numerical Methods

In this section, we will introduce the DG algorithm for the VM system. We consider an infinite, homogeneous plasma, where all boundary conditions in are periodic, and is assumed to be compactly supported in . This assumption is consistent with the fact that the solution of the VM system is expected to decay at infinity in -space, preserving integrability and its kinetic energy.

Without loss of generality, we assume and , where the velocity space domain is chosen large enough so that at and near the phase space boundaries. We take in the following sections, although the method and its analysis can be extended directly to the cases when and take any values from .

In our analysis, the assumption that remain compactly support in , given that it is initially so, is an open question in the general setting. The answer to this question is important for proving the existence of a globally defined classical solution, and its failure could indicate the formation of shock-like solutions of the VM system. Whether or not the three-dimensional VM system is globally well-posed as a Cauchy problem is a major open problem. The limited results of global existence without uniqueness of weak solutions and well-posedness and regularity of solutions assuming either some symmetry or near neutrality constitute the present extent of knowledge [29, 30, 25, 21, 26, 28, 27].

### 2.1 Notations

Throughout the paper, standard notations will be used for the Sobolev spaces. Given a bounded domain (with , or ) and any nonnegative integer , denotes the -Sobolev space of order with the standard Sobolev norm , and denotes the -Sobolev space of order with the standard Sobolev norm and the semi-norm . When , we also use and .

Let and be partitions of and , respectively, with and being (rotated) Cartesian elements or simplices; then defines a partition of . Let be the set of the edges of and be the set of the edges of ; then the edges of will be . Here we take into account the periodic boundary condition in the -direction when defining and . Furthermore, with and being the set of interior and boundary edges of , respectively. In addition, we denote the mesh size of as , where with , with , and for . When the mesh is refined, we assume both and are uniformly bounded from above by a positive constant . Here and . It is further assumed that is shape-regular with or . That is, if denotes the diameter of the largest sphere included in , there is

 hK⋆ρK⋆≤σ⋆,∀K⋆∈T⋆h

for a positive constant independent of .

Next we define the discrete spaces

 Gkh ={g∈L2(Ω):g|K=Kx×Kξ∈Pk(Kx×Kξ),∀Kx∈Txh,∀Kξ∈Tξh} , (2.2a) ={g∈L2(Ω):g|K∈Pk(K),∀K∈Th} , Urh ={U∈[L2(Ωx)]dx:U|Kx∈[Pr(Kx)]dx,∀Kx∈Txh} , (2.2b)

where denotes the set of polynomials of total degree at most on , and and are nonnegative integers. Note the space , which we use to approximate , is called P-type, and it can be replaced by the tensor product of P-type spaces in and ,

 {g∈L2(Ω):g|K=Kx×Kξ∈Pk(Kx)×Pk(Kξ),∀Kx∈Txh,∀Kξ∈Tξh} , (2.3)

or by the tensor product space in each variable, which is called Q-type

 {g∈L2(Ω):g|K=Kx×Kξ∈Qk(Kx)×Qk(Kξ),∀Kx∈Txh,∀Kξ∈Tξh} . (2.4)

Here denotes the set of polynomials of degree at most in each variable on . The numerical methods formulated in this paper, as well as the conservation, stability, and error estimates, hold when any of the spaces above is used to approximate . In our simulations of Section 5, we use the P-type of (2.2a) as it is the smallest and therefore renders the most cost efficient algorithm. In fact, the ratios of these three spaces defined in (2.2a), (2.3) and (2.4) are with .

For piecewise functions defined with respect to or , we further introduce the jumps and averages as follows. For any edge , with as the outward unit normal to , , and , the jumps across are defined as

 [g]x=g+n+x+g−n−x,[U]x=U+⋅n+x+U−⋅n−x,[U]τ=U+×n+x+U−×n−x

and the averages are

 {g}x=12(g++g−),{U}x=12(U++U−).

By replacing the subscript with , one can define , , , and for an interior edge of in . For a boundary edge with being the outward unit normal, we use

 [g]ξ=gnξ,{g}ξ=12g,{U}ξ=12U . (2.5)

This is consistent with the fact that the exact solution is compactly supported in .

For convenience, we introduce some shorthand notations, where again is or . In addition, with and . There are several equalities that will be used later, which can be easily verified using the definitions of averages and jumps.

 12[g2]⋆={g}⋆[g]⋆,with⋆=xorξ , (2.6a) [U×V]x+{V}x⋅[U]τ−{U}x⋅[V]τ=0 , (2.6b) [U×V]x+V+⋅[U]τ−U−⋅[V]τ=0,[U×V]x+V−⋅[U]τ−U+⋅[V]τ=0 . (2.6c)

We end this subsection by summarizing some standard approximation properties of the above discrete spaces, as well as some inverse inequalities [14]. For any nonnegative integer , let be the projection onto , and be the projection onto , then

###### Lemma 2.1 (Approximation properties)

There exists a constant , such that for any and , the following hold:

 ||g−Πmg||0,K+h1/2K||g−Πmg||0,∂K ≤Chm+1K||g||m+1,K,∀K∈Th , ||U−\boldmath{Π}mxU||0,Kx+h1/2Kx||U−\boldmath{Π}mxU||0,∂Kx ≤Chm+1Kx||U||m+1,Kx,∀Kx∈Txh , ||U−\boldmath{Π}mxU||0,∞,Kx ≤Chm+1Kx||U||m+1,∞,Kx,∀Kx∈Txh ,

where the constant is independent of the mesh sizes and , but depends on and the shape regularity parameters and of the mesh.

###### Lemma 2.2 (Inverse inequality)

There exists a constant , such that for any or with , and for any , the following hold:

 ||∇xg||0,K≤Ch−1Kx||g||0,K,||∇ξg||0,K≤Ch−1Kξ||g||0,K,
 ||U||0,∞,Kx≤Ch−dx/2Kx||U||0,Kx,||U||0,∂Kx≤Ch−1/2Kx||U||0,Kx ,

where the constant is independent of the mesh sizes , , but depends on and the shape regularity parameters and of the mesh.

### 2.2 The Semi-Discrete DG Methods

On the PDE level, the two equations in (1.1c) involving the divergence of the magnetic and electric fields can be derived from the remaining part of the VM system; therefore, the numerical methods proposed in this section are formulated for the VM system without (1.1c). We want to stress that even though in principle the initial satisfaction of these constraints is sufficient for their satisfaction for all time, in certain circumstance one may need to consider explicitly such divergence conditions in order to produce physically relevant numerical simulations [46, ?].

Given , the semi-discrete DG methods for the VM system are defined by the following procedure: for any , look for , , such that for any , ,

 ∫K∂tfhgdxdξ −∫Kfhξ⋅∇xgdxdξ−∫Kfh(Eh+ξ×Bh)⋅∇ξgdxdξ +∫Kξ∫∂Kxˆfhξ⋅nxgdsxdξ+∫Kx∫∂Kξˆ(fh(Eh+ξ×Bh)⋅nξ)gdsξdx=0 , (2.7a) ∫Kx∂tEh⋅Udx =∫KxBh⋅∇×Udx+∫∂Kxˆnx×Bh⋅Udsx−∫KxJh⋅Udx , (2.7b) ∫Kx∂tBh⋅Vdx =−∫KxEh⋅∇×Vdx−∫∂Kxˆnx×Eh⋅Vdsx , (2.7c)

with

 Jh(x,t)=∫Tξhfh(x,ξ,t)ξdξ . (2.8)

Here and are outward unit normals of and , respectively. All ‘hat’ functions are numerical fluxes that are determined by upwinding, i.e.,

 ˆfhξ⋅nx: =˜fhξ⋅nx=({fhξ}x+|ξ⋅nx|2[fh]x)⋅nx , (2.9a) ˆfh(Eh+ξ×Bh)⋅nξ: =˜fh(Eh+ξ×Bh)⋅nξ =({fh(Eh+ξ×Bh)}ξ+|(Eh+ξ×Bh)⋅nξ|2[fh]ξ)⋅nξ , (2.9b) ˆnx×Eh: =nx×˜Eh=nx×({Eh}x+12[Bh]τ) , (2.9c) ˆnx×Bh: =nx×˜Bh=nx×({Bh}x−12[Eh]τ) , (2.9d)

where these relations define the meaning of ‘tilde’.

For the Maxwell part, we also consider two other numerical fluxes: central flux and alternating flux, which are defined by

 Central flux:˜Eh={Eh},˜Bh={Bh} , (2.10a) Alternating flux:˜Eh=E+h,˜Bh=B−h,or˜Eh=E−h,˜Bh=B+h . (2.10b)

Upon summing up (2.7a) with respect to and similarly summing (2.7b) and (2.7c) with respect to , the numerical method becomes the following: look for , , such that

 ah(fh,Eh,Bh;g) =0 , (2.11a) bh(Eh,Bh;U,V) =lh(Jh;U) , (2.11b)

for any , , where

 ah(fh,Eh,Bh;g)= ah,1(fh;g)+ah,2(fh,Eh,Bh;g) ,lh(Jh;U)=−∫TxhJh⋅Udx bh(Eh,Bh;U,V)= ∫Txh∂tEh⋅Udx−∫TxhBh⋅∇×Udx−∫Ex˜Bh⋅[U]τdsx +∫Txh∂tBh⋅Vdx+∫TxhEh⋅∇×Vdx+∫Ex˜Eh⋅[V]τdsx ,

and

 ah,1(fh;g) =∫Th∂tfhgdxdξ−∫Thfhξ⋅∇xgdxdξ+∫Tξh∫Ex˜fhξ⋅[g]xdsxdξ , ah,2(fh,Eh,Bh;g) =−∫Thfh(Eh+ξ×Bh)⋅∇ξgdxdξ+∫Txh∫Eξ˜fh(Eh+ξ×Bh)⋅[g]ξdsξdx .

Note, is linear with respect to and , yet it is in general nonlinear with respect to and due to (2.9b). Recall, the exact solution has compact support in ; therefore, the numerical fluxes of (2.9a)-(2.9d) and (2.10a) and (2.10b) are consistent and, consequently, so is the proposed numerical method. That is, the exact solution satisfies

 ah(f,E,B;g) =0,∀g∈Gkh , bh(E,B;U,V) =lh(J;U),∀U,V∈Urh .

### 2.3 Temporal Discretizations

We use total variation diminishing (TVD) high-order Runge-Kutta methods to solve the method of lines ODE resulting from the semi-discrete DG scheme, . Such time stepping methods are convex combinations of the Euler forward time discretization. The commonly used third-order TVD Runge-Kutta method is given by

 G(1)h = Gnh+△tR(Gnh) G(2)h = 34Gnh+14G(1)h+14△tR(G(1)h) Gn+1h = 13Gnh+23G(2)h+23△tR(G(2)h), (2.12)

where represents a numerical approximation of the solution at discrete time . A detailed description of the TVD Runge-Kutta method can be found in [53]; see also [31] and [32] for strong-stability-perserving methods.

## 3 Conservation and Stability

In this section, we will establish conservation and stability properties of the semi-discrete DG methods. In particular, we prove that subject to boundary effects, the total charge (mass) is always conserved. As for the total energy of the system, conservation depends on the choice of numerical fluxes for the Maxwell equations. We also show that is stable, which facilitates the error analysis of Section 4.

###### Lemma 3.1 (Mass conservation)

The numerical solution with satisfies

 ddt∫Thfhdxdξ+Θh,1(t)=0 , (3.13)

where

 Θh,1(t)=∫Txh∫Ebξfhmax((Eh+ξ×Bh)⋅nξ,0)dsξdx .

Equivalently, with , for any , the following holds:

 ∫Txhρh(x,T)dx+∫T0Θh,1(t)dt=∫Txhρh(x,0)dx . (3.14)
{proof}

Let . Note that , for any , is continuous and . Taking this as the test function in (2.11a), one has

 ddt∫Thfhdxdξ+∫Txh∫Ebξ˜fh(Eh+ξ×Bh)⋅[g]ξdsξdx=0 .

With the numerical flux of (2.9b) and the average and jump across of (2.5), the second term above becomes

 ∫Txh∫Ebξ˜fh(Eh+ξ×Bh)⋅nξdsξdx (3.15) = ∫Txh∫Ebξfh2((Eh+ξ×Bh)⋅nξ+|(Eh+ξ×Bh)⋅nξ|)dsξdx=Θh,1(t) , (3.16)

and this gives (3.13). Integrating in time from to gives (3.14).

###### Lemma 3.2 (Energy conservation 1)

For , , the numerical solution , with the upwind numerical fluxes (2.9a)-(2.9d) satisfies

 ddt(∫Thfh|ξ|2dxdξ+∫Txh(|Eh|2+|Bh|2)dx)+Θh,2(t)+Θh,3(t)=0 , (3.17)

with

 Θh,2(t)=∫Ex(|[Eh]τ|2+|[Bh]τ|2)dsx ,Θh,3(t)=∫Txh∫Ebξfh|ξ|2max((Eh+ξ×Bh)⋅nξ,0)dsξdx .
{proof}

Step 1: Let . Note that for and it is continuous. In addition, , , and for any function . Taking this as the test function in (2.11a), one has

 ddt∫Thfh|ξ|2dxdξ=2∫ThfhEh⋅ξdxdξ−∫Txh∫Ebξ˜fh(Eh+ξ×Bh)⋅[|ξ|2]ξdsξdx =2∫TxhEh⋅(∫Tξhfhξdξ)dx−∫Txh∫Ebξ(12(Eh+ξ×Bh)fh+|(Eh+ξ×Bh)⋅nξ|2fhnξ)⋅(|ξ|2nξ)dsξdx =2∫TxhEh⋅Jhdx−∫Txh∫Ebξfh2((Eh+ξ×Bh)⋅nξ+|(Eh+ξ×Bh)⋅nξ|)|ξ|2d