High Order Maximum Principle Preserving SemiLagrangian Finite Difference WENO schemes for the Vlasov Equation
Tao Xiong ^{1}^{1}1Department of Mathematics, University of Houston, Houston, 77204. Email: txiong@math.uh.edu. JingMei Qiu ^{2}^{2}2Department of Mathematics, University of Houston, Houston, 77204. Email: jingqiu@math.uh.edu. The first and second authors are supported by Air Force Office of Scientific Computing YIP grant FA9550120318, NSF grant DMS1217008 and University of Houston. Zhengfu Xu ^{3}^{3}3Department of Mathematical Science, Michigan Technological University, Houghton, 49931. Email: zhengfux@mtu.edu. Supported by NSF grant DMS1316662. Andrew Christlieb ^{4}^{4}4 Department of Mathematics, Michigan State University, East Lansing, MI, 48824. Email: andrewch@math.msu.edu
Abstract
In this paper, we propose the parametrized maximum principle preserving (MPP) flux limiter, originally developed in [Z. Xu, Math. Comp., (2013), in press], to the semiLagrangian finite difference weighted essentially nonoscillatory scheme for solving the Vlasov equation. The MPP flux limiter is proved to maintain up to fourth order accuracy for the semiLagrangian finite difference scheme without any time step restriction. Numerical studies on the VlasovPoisson system demonstrate the performance of the proposed method and its ability in preserving the positivity of the probability distribution function while maintaining the high order accuracy.
Keywords: SemiLagrangian method; Finite difference WENO scheme; Maximum principle preserving; Parametrized flux limiter; Vlasov equation
1 Introduction
In this paper, we will consider the VlasovPoisson (VP) system
(1.1)  
(1.2) 
on the domain , where and . It is an important system for modelling the collisionless plasma. The Vlasov equation (1.1) is a kinetic equation that describes the time evolution of the probability distribution function (PDF) of finding an electron at position with velocity at time . is the electric field and is a selfconsistent electrostatic potential function described by the Poisson equation (1.2). The probability distribution function is coupled with the long range fields via the charge density , with uniformly distributed infinitely massive ions in the background. The total charge neutrality condition is imposed. The Vlasov equation (1.1) in the hyperbolic form enjoys the maximum principle preserving (MPP) and positivity preserving (PP) properties. That is, let
the solution at later time is still within the interval and stays positive.
The VP system has been extensively studied numerically. Popular numerical methods, besides the semiLagrangian (SL) approaches which will be reviewed in the next paragraph, include particleincell (PIC) methods [12, 5, 20], Lagrangian particle methods [4, 13], weighted essentially nonoscillatory (WENO) coupled with Fourier collocation [42], FourierFourier spectral methods [23, 22], finite volume methods [14, 16], continuous finite element methods [39, 38], RungeKutta discontinuous Galerkin methods [19, 3, 11, 9, 8], and methods in references therein.
In this paper, we focus on the SL method, which has been proposed and applied for a wide range of applications, e.g. in atmospheric modeling and simulations [32, 17], in capturing the moving interface via solving level set equations [33], in fluid dynamics [36] and in kinetic simulations [15, 34]. Compared to the Eulerian approach, the SL approach also has a fixed set of numerical meshes, but in each time step evolution, the information propagates along the characteristics in the Lagrangian fashion. The SL method can be designed as accurate as the Eulerian approach. In addition, the SL method is free of time step restriction by utilizing the characteristic method in the temporal direction. The design of SL approach requires the tracking of characteristics forward or backward in time, which could be challenging for nonlinear problems. To avoid such difficulty in solving the VP system, the dimensional splitting approach was proposed in [7]. Along this line, SL methods were developed in various settings. For example, in the finite difference framework, different interpolation strategies, such as the cubic spline interpolation [31], the cubic interpolated propagation [25], the high order WENO interpolation in a nonconservative form [6] and in a conservative form [26, 27, 28] are proposed. There have also been many work in designing the SL methods in the finite volume framework for the VP system [16] and for the guiding center Vlasov model [10], in the finite element discontinuous Galerkin (DG) framework [30, 29] and with a hybrid finite differencefinite element approach [18].
Our main interest in this paper is to develop a MPP, and in particular PP, SL finite difference WENO scheme for the VP system (1.1). The main challenge of designing MPP (and/or PP) schemes within the SL WENO framework is to maintain the designed high order of accuracy from the conservative WENO approximation. Meanwhile, it is desired that no additional restrictive CFL constraint will be introduced. There have been some research efforts on designing MPP and PP high order SL schemes. For example, in [16, 10], PP property is preserved in the finite volume framework. However, the PP preservation of the PDF is accompanied with sacrificing high order spatial accuracy. Recently, in [30, 29], the SL discontinuous Galerkin (DG) methods to solve the VP system is coupled with the MPP limiters, that were originally proposed by Zhang and Shu [40]. In these approaches, the limiters are applied to the reconstructed polynomials (for finite volume) or the representing polynomials (for DG). In general, the MPP (or PP) property together with the maintenance of high order accuracy is much more challenging to achieve in the finite difference framework than in the finite volume and finite element DG framework via limiting polynomials, see [41].
We propose to generalize the recently developed parametrized MPP flux limiter [37] to a conservative SL finite difference WENO method solving the VP system. The original parametrized flux limiter for 1D scalar conservation laws was later extended to the twodimensional case [24]. Xiong et al. [35] proposed to apply the parametrized MPP flux limiter to the final RungeKutta (RK) stage only, with significantly improved time step restriction for maintenance of high order accuracy, leading to much reduced computational cost. It has also been proved in [35] that the parametrized MPP flux limiter can maintain up to third order accuracy both in space and in time for nonlinear scalar conservation laws.
In order to apply the parametrized MPP flux limiter in [37] for solving the VP system, we start with the dimensional splitting approach. The parametrized MPP flux limiter is proposed for the conservative high order SL finite difference WENO scheme [27, 26]. We mimic the proof in [35] to prove that the parametrized MPP flux limiter for the SL finite difference scheme solving linear advection equations can maintain up to fourth order accuracy without any time step constraint. We also apply the parametrized flux limiter proposed in [35] to the RK finite difference WENO approximation of the VP system without operator splitting. Through numerical studies on weak and strong Landau damping, two stream instabilities, and KEEN waves, we show that both methods perform very well with the designed MPP properties, while maintaining high order accuracy and mass conservation.
The rest of the paper is organized as follows. In Section 2, the SL finite difference WENO scheme is reviewed and the parametrized MPP flux limiting procedure is proposed for the SL WENO scheme. We also prove that the parametrized MPP flux limiter maintains up to fourth order accuracy without additional time step restriction. Numerical studies of the scheme are presented in Section 3. Conclusions are made in Section 4.
2 The parametrized MPP flux limiter for the SL WENO scheme
In this section, we will briefly describe the SL finite difference WENO scheme for solving Vlasov equation (1.1) in one dimension (both physical and phase spaces). We adopt the method in [28], which is in a conservative flux difference form and is suitable to be coupled with the newly developed parametrized MPP flux limiter [37].
For the onedimensional problem, periodic boundary conditions are imposed in xdirection on a finite domain and in vdirection with a cutoff domain with chosen to be large enough to guarantee for . The domain is discretized by the the computational grid
Let and to be the middle point of each cell. The uniform mesh sizes in each direction are and .
Following [7], the Vlasov equation is dimensionally split to the following form
(2.1)  
(2.2) 
Using the second order Strang splitting strategy, the numerical solution is updated from time level to time level by solving equation (2.1) for half a time step, then solving equation (2.2) by a full time step, followed by solving equation (2.1) for another half a time step. Note that each split equation (2.1) or (2.2) still preserves the MPP (and/or PP) property.
In the following, we will take the prototype 1D linear advection equation
(2.3) 
with constant , to present the SL finite difference scheme with the parametrized MPP flux limiters.
2.1 Review of SL finite difference WENO scheme
The SL finite difference WENO scheme proposed in [28] is based on integrating equation (2.3) in time over ,
(2.4) 
where
(2.5) 
By introducing a sliding average function
(2.6) 
with
(2.7) 
the evaluation of equation (2.4) at the grid point can be written in a conservative form
(2.8) 
where is called the flux function. A SL finite difference WENO reconstruction is used to approximate the numerical flux function based on its cell averages
(2.9) 
with
where is the point tracing from the grid point along characteristics back to the time level . The last equality above is essential for the SL scheme, and is obtained via following characteristics. can be reconstructed from for each . In summary, the SL finite difference WENO scheme procedure in evolving equation (2.3) from to is as follows:

At each of the grid points at time level , say , trace the characteristic back to time level at .

Reconstruct from . We use to denote this reconstruction procedure
(2.10) to approximate , where indicates the stencil used in the reconstruction. indicates the reconstruction of .

Reconstruct from . We use to denote this reconstruction procedure
(2.11) to approximate . Here indicates the stencil used in the reconstruction.

Update the solution by
(2.12) with numerical fluxes computed in the previous step.
When the reconstruction stencils in and above only involve one neighboring point value of the solution, then the scheme reduces to a first order monotone scheme when the time step is within CFL restriction. We let denote the first order flux. The proposed SL finite difference scheme can be designed to be of high order accuracy by including more points in the stencil for (the composition of and ), to reconstruct the numerical flux
(2.13) 
where indicates the stencil used in the reconstruction process. The WENO mechanism can be introduced in the reconstruction procedures in order to realize a stable and nonoscillatory capture of fine scale structures. In the Appendix, we provide formulas to obtain the high order fluxes in (2.13) for the fifth order SL finite difference WENO scheme. For more details, we refer to [28].
For the case of large time step , if , is no longer inside , let to be the index such that and , we have
(2.14)  
(2.15) 
where is reconstructed in the same fashion as (2.13), but replacing by . Similarly, if , let to be the index such that and , we have
(2.16)  
(2.17) 
It is numerically demonstrated in [28] that the proposed high order SL WENO method works very well in Vlasov simulations with extra large time step evolution.
2.2 Parametrized MPP flux limiters
In this subsection, we propose a parametrized MPP flux limiter, as proposed in [37, 35], for the high order SL finite difference WENO scheme (2.12).
For simplicity, Let and as the minimum and maximum values of the initial condition. It has been known that the numerical solutions updated by (2.12) with the first order monotone flux satisfy the maximum principle. The MPP flux limiters are designed as modifying the high order numerical flux towards the first order monotone flux in the following way,
(2.18) 
where is the parameter to be designed to take advantage of the first order monotone flux in the MPP property and to take advantage of the high order flux in the high order accuracy.
Below is a detailed procedure of designing , in order to guarantee the MPP property of the numerical solutions, yet to choose ’s to be as close to as possible for high order accuracy. For each in limiting the numerical flux as in equation (2.18), we are looking for the upper bounds and , such that, for all
(2.19) 
the updated numerical solutions and by the SL WENO scheme (2.12) with the modified numerical fluxes (2.18) are within . Let
the MPP property of the first order monotone flux guarantees
(2.20) 
To ensure with as in equation (2.18), it is sufficient to have
(2.21)  
(2.22) 
The linear decoupling of , subject to the constraints (2.21) and (2.22), is achieved via a casebycase discussion based on the sign of
as outlined below.
Notice that the range of (2.19) is determined by the need to ensure both the upper bound (2.21) and the lower bound (2.22) of numerical solutions in both cell and . Thus the locally defined limiting parameter is given as
(2.23) 
with and . The modified flux in equation (2.18) with the designed above ensures the maximum principle. Such modified flux is consistent and preserves the maximum principle by its design, since it is a convex combination () of a high order flux and the first order flux . The mass conservation property is preserved, due to the flux difference form of the scheme.
Remark 2.1.
In [35], it was proved that the MPP flux limiter for the RK finite difference scheme can maintain up to third order accuracy both in space and in time for the nonlinear scalar equation . In Theorem 2.2 below, we prove the proposed MPP limiter for the SL finite difference scheme can maintain up to fourth order accuracy for solving the linear advection equation (2.3) without any time step restriction. The generalization of the proof to maintain up to fifth order accuracy would be much more technical and will be investigated in our future work.
Theorem 2.2.
Proof.
Without loss of generality, we consider with . The case of can be reduced to by shifting the numerical solution with several whole grid points. Without specifying, in the following, we use instead of and use instead of . Since the difference between and is of high order due to assumption (2.24), we use and interchangeably when such high order difference allows. The high order flux is taken to be (2.13) with a 4th order finite difference reconstruction, the first order monotone flux is .
We mimic the proof in [35] and only consider the most challenging case: case (b) for the maximum value part. The proof for other cases would be similar to that in [35].
We consider case (b) when
(2.26) 
To prove (2.25), it suffices to prove
(2.27) 
if . For the SL method, we have
(2.28) 
where is defined by
(2.29) 
follows the same characteristics as the linear advection equation, hence
(2.30) 
Let . We approximate by a 4th order reconstructed polynomial from the solution based on the stencil , with which
(2.31)  
We first consider the case when the maximum point at , with , and . is the dependent domain for the exact solution when . We perform Taylor expansions around up to th order and obtain
Since , we can write to be
(2.32) 
where
and .
We consider a quantity denoted by , which can be written as follows
(2.33)  
with parameters and to be determined. If , .
Since , if we can find , such that
(2.34) 
we would have
(2.35) 
when (2.32) is compared with (2.33) under the assumption (2.26), thus proving (2.27) in this case. In the following, to determine the parameters and satisfying (2.34), we first need , that is
(2.36) 
Using Mathematica, we have
(2.37) 
and
(2.38)  
(2.39) 
Thus, we can choose . Let . Below, we verify that is bounded.