1 Introduction

A WENO-based Method of Line Transpose Approach for Vlasov Simulations

Andrew Christlieb111 Department of Mathematics and Electrical Engineering, Michigan State University, East Lansing, MI, 48824. E-mail: andrewch@math.msu.edu , Wei Guo 222 Department of Mathematics, Michigan State University, East Lansing, MI, 48824. E-mail: wguo@math.msu.edu, Yan Jiang333Department of Mathematics, Michigan State University, East Lansing, MI, 48824. E-mail: jiangyan@math.msu.edu

Abstract.

In this paper, a high order implicit Method of Line Transpose (MOL) method based on a weighted essentially non-oscillatory (WENO) methodology is developed for one-dimensional linear transport equations and further applied to the Vlasov-Poisson (VP) simulations via dimensional splitting. In the MOL framework, the time variable is first discretized by a diagonally implicit strong-stability-preserving Runge-Kutta method, resulting in a boundary value problem (BVP) at the discrete time levels. Then an integral formulation coupled with a high order WENO methodology is employed to solve the BVP. As a result, the proposed scheme is high order accurate in both space and time and free of oscillations even though the solution is discontinuous or has sharp gradients. Moreover, the scheme is able to take larger time step evolution compared with an explicit MOL WENO scheme with the same order of accuracy. The desired positivity-preserving (PP) property of the scheme is further attained by incorporating a newly proposed high order PP limiter. We perform numerical experiments on several benchmarks including linear advection, solid body rotation problem; and on the Landau damping, two-stream instabilities, bump-on-tail, and plasma sheath by solving the VP system. The efficacy and efficiency of the proposed scheme is numerically verified.

Key Words: Method of Line Transpose, implicit Runge-Kutta methods, weighted essentially non-oscillatory methodology, high order accuracy, Vlasov-Poisson, demensional splitting, positivity-preserving

## 1 Introduction

Transport modeling plays an important role in a wide range of applications, including fluid dynamics, plasma physics, atmospheric sciences, among many others. The main focus of this paper is to develop of efficient high order transport solvers for the Vlasov simulations. As a high-dimensional kinetic transport equation, the Vlasov-Poisson (VP) system is considered a fundamental model in plasma physics which describes the dynamics of charged particles due to the self-consistent electric force:

 ft+v⋅∇xf+E(x,t)⋅∇vf=0,x×v∈Ωx×Ωv (1.1a) E(x,t)=−∇xϕ(x,t),  −Δxϕ(x,t)=−1+ρ(x,t) (1.1b)

where describes the probability of finding a particle with velocity at position at time , is the electrostatic field, is the self-consistent electrostatic potential, and is the electron charge density and the 1 represents the uniformly distributed infinitely massive ions in the background.

The VP simulation is challenging. Besides the famous curse of dimensionality, the underlying system may develop filamentation solution structures, since the long-range forces drive the solution out of thermo-equilibrium. Hence, it is of importance that the used numerical scheme is able to effectively capture filamentation structures without producing spurious oscillations. Further, the VP system conserves several physical quantities including the total mass, momentum, and energy, which requires that the numerical scheme converses those invariants on the discrete level. In addition, as a probability density function in phase space, the unknown function always remains non-negative. In practice, the positivity preserving property is another highly desired property for the design of Vlasov solver, since negative solutions may trigger unphysical phenomenon in plasmas.

There exist a large amount of successful VP solvers in the literature, such as the prominent particle-in-cell (PIC) method [4, 23, 41] , which evolves the solution by following the trajectories of some sampled macro-particles, Eulerian type approaches [18, 50, 12, 11], for which the state variable is advanced according to the underlying equations on a fixed numerical grid, and semi-Lagrangian approach [9, 29, 2, 39, 19, 3, 17, 16, 34, 30, 15, 22], which updates the solution by tracking the characteristics on a fixed mesh grid. However, these popular schemes still have limitations. For example, the PIC method is known to suffer from low order sampling noise; Eulerian approaches are subject to stringent time step restriction associated with explicit time stepping; general boundary conditions are difficult to handle in a semi-Lagrangian setting. Also note that many of these existing methods are in the Method of Line (MOL) framework, meaning that the spatial variable is first discretized, then the numerical solution is updated in time by coupling a suitable time integrator. In this paper, we consider an alternative approach to advance the solution, which is called the Method of Line Transpose (MOL) or Rothes’s method in the literature [35, 36]. As its name implies, such an approach is defined in an orthogonal fashion of the MOL schemes, i.e., the discretization is first carried out for temporal variable, resulting in a boundary value problem (BVP) at the discrete time levels. Then a preferred BVP solver can be further applied to advance the solution. For example, in [6], an implicit integral formulation based on the MOL framework for one-dimensional wave equations is developed and analyzed. It is proven that the scheme is A-stable, hence allowing for extra large time step evolution. In [8], the computational cost of the scheme is reduced from to , where is the number of discrete mesh points, by taking advantage of an important property of the analytical solution of the wave equations. More recently, a successive convolution technique [5] is developed to enhance the temporal accuracy to arbitrary high order. A notable advantage of such a scheme is that, even though it is implicit in time, we do not need to explicitly solve a linear system, while the BVP is inverted analytically in an integral formulation. The method is further extended to multi-dimensional wave equations by utilizing the alternating direction implicit strategy at the expense of introducing splitting errors [6]. However, the simplicity of the scheme is preserved. The recent applications of this MOL approach include the heat equation, the Allen-Cahn equation [7], and Maxwell’s equations [10].

In this paper, we consider formulating a numerical scheme suited for the Vlasov equation in the MOL framework in order to take advantage of its good properties such as the simplicity and the ability to take large time steps for evolution. In order to accommodate the MOT approach, the nonlinear multi-dimensional Vlasov equation is first split into several lower-dimensional linear transport equations. Such an idea has been widely used in semi-Lagrangian Vlasov simulations [9]. As we mentioned earlier, the extra difficulty of VP simulations lies in the distinctive filamentation solution structures, which requires that the employed numerical scheme is able to effectively capture sharp gradient of the solution without producing spurious oscillations even on an under-resolved mesh. Note that we cannot directly apply the existing MOL schemes to the VP system, since they are all linear schemes and hence exhibit oscillations when the solution is discontinuous or has a sharp gradient. In order to control undesired oscillations, we propose to further couple two numerical ingredients. Firstly, in the MOL framework, the time variable is discretized via an implicit high order strong-stability-preserving (SSP) Runge-Kutta (RK) method [21]. Secondly, a robust weighted essential non-oscillatory (WENO) methodology is incorporated when inverting the BVP. We expect that the resulting scheme is efficient, high order accurate, robust, and able to take large time step evolution (larger than the explicit counterpart). Furthermore, a high order positivity-preserving (PP) limiter suited for the proposed MOL scheme is developed to ensure such a high desired property for the numerical solution.

We also would like to point out that the proposed MOL WENO scheme shares some similarity with the semi-Lagrangian WENO scheme developed in [31, 32], in the sense that both schemes are high order accurate and able to effectively capture the filamentation structures with large time step. On the other hand, compared with semi-Lagrangian WENO schemes, the proposed method is able to conveniently handle general boundary conditions, which is not trivial for semi-Lagrangian approaches.

The rest of the paper is organized as follows. In Section 2, we formulate the numerical scheme for one-dimensional linear transport equations by coupling several main numerical ingredients such as the MOL approach, the WENO methodology and the high order PP limiter. In Section 3, the method is extended to two-dimensional transport equations via dimensional splitting. We present several benchmark tests in Section 4 to verify the performance of the proposed scheme, including linear transport, rigid body rotation, as well as Landau damping, two-stream instabilities, bump-on-tail and plasma sheath for VP simulations. We conclude the paper in Section 5.

## 2 Formulation of the Scheme

 ut+cux=0,x∈[a,b], (2.2)

where is constant and denotes the wave propagation speed. The boundary condition is assumed to be periodic, Dirichlet boundary condition

 u(a,t)=g1(t), for c>0,   or    u(b,t)=g2(t), for c<0,

or Neumann boundary condition

 ux(a,t)=h1(t), for c>0,   or    ux(b,t)=h2(t), for c<0.

### 2.1 The method of line transpose approach

As discussed in the introduction, in a MOL framework, we first discretize the time variable while leaving the spatial variable continuous. For example, we can apply the first order backward Euler method to (2.2), and obtain

 un+1−unΔt+cun+1x=0,

i.e.,

 un+1x+sgn(c)⋅αun+1=sgn(c)⋅αun, (2.3)

where denotes the time step and . Note that (2.3) is an initial value problem (IVP), for which we can obtain an explicit representation for solution at time level . If we assume , can be represented by

 un+1(x)=IL[un,α](x)+An+1e−α(x−a), (2.4)

where

 IL[v,α](x)≐α∫xae−α(x−y)v(y)dy, (2.5)

and is a constant and should be determined by the boundary condition. The superscript indicates that the characteristics traverse from the left to the right. If , then the characteristics traverse from the right to the left. Due to the symmetry, we can advance the solution as

 un+1(x)=IR[un,α](x)+Bn+1e−α(b−x), (2.6)

when the first order backward Euler scheme is used. Here,

 IR[v,α](x)≐α∫bxe−α(y−x)v(y)dy, (2.7)

and again the constant is determined by the boundary condition.

We need to further discretize the spatial variable based on the semi-discrete schemes (2.4) and (2.6). Following the idea in [6, 8], the domain is discretized by uniformly distributed grid points

 a=x0

and denotes the mesh size. The solution is accordingly discretized at points , and we use to denote the numerical solution at location at time level . The main part of the algorithm is to evaluate the convolution integral in (2.4), (2.6) or in (2.5), (2.7) based on the point values .

We first consider the case of . Note that,

 IL[v,α](x)=IL[v,α](x−δx)e−αδx+JL[v,α,δx](x), (2.8)

where

 JL[v,α,δx](x)≐α∫xx−δxv(y)e−α(x−y)dy, (2.9)

and is a positive constant. Therefore, in order to obtain denoted by , we only need to approximate the increment denoted by , and obtain via a recursive relation, i.e.,

 ILi=ILi−1e−αΔx+JLi,i=1,⋯,M,IL0=0. (2.10)

For the case of , again thanks to the symmetry, we have the following identity

 IR[v,α](x)=IR[v,α](x+δx)e−αδx+JR[v,α,δx](x), (2.11)

where

 JR[v,α,δx](x)≐α∫x+δxxv(y)e−α(y−x)dy, (2.12)

and hence a similar recursive relation can be obtained

 IRi=IRi+1e−αΔx+JRi,i=0,⋯,M−1,IRM=0, (2.13)

where and denote the approximate values of and , respectively.

After computing the contributions from the left or right characteristic, we further need to find the global coefficient in (2.4) or in (2.6), which is used to enforce the boundary condition. Since , for a Dirichlet boundary condition, the global coefficient can be obtained by

 An+1=g1(tn+1), or  Bn+1=g2(tn+1). (2.14)

Similarly, for a Neumann boundary condition, because and , we can get the coefficients by

 An+1=v(a)−1αh1(tn+1),  or  Bn+1=v(b)+1αh2(tn+1). (2.15)

For the periodic boundary condition, the coefficient can be determined by using the mass conservation property of the solution, i.e,

 ∫bau(x,t)dx=∫bau(x,0)dx.

On the discrete level, we enforce

 M−1∑i=0uni=M−1∑i=0un+1i=M−1∑i=0ILi+An+1M−1∑i=0e−iαΔx,   un+1M=un+10,   for c>0,

or

 M∑i=1uni=M∑i=1un+1i=M∑i=1IRi+Bn+1M∑i=1e−(M−i)αΔx,   un+10=un+1M,   for c<0.

Then the coefficient is obtained

 (2.16)

or

 Bn+1=(M∑i=1uni−M∑i=1IRi)/M∑i=1e−(M−i)αΔx. (2.17)

### 2.2 WENO methodology

The existing MOL schemes mainly use linear interpolation methods to compute increment and such as [6], which work well for smooth problems. However, the schemes may generate spurious oscillations when the solution has discontinuities because of the famous Gibbs phenomenon. In order to control such undesired oscillations, we will incorporate a WENO methodology in the MOL framework to evaluate and below.

The first WENO scheme was introduced by Liu et al. [27], as a third order accurate finite volume scheme to solve hyperbolic conservation laws in one dimension. Later, in [24], a general framework for arbitrary order accurate finite difference WENO schemes was designed, which is more efficient for multi-dimensional calculations. The schems have the advantage of attaining high order accuracy in smooth regions of the solution while maintaining sharp and essentially non-oscillatory transitions of discontinuities. More information about WENO can be found in [38].

In our work, a WENO integration procedure is used, based on the idea proposed in [13, 14, 28]. Below, a -th order approximation of is provided as an example, where is an integer. The corresponding coefficients of the third order and fifth order schemes are given in the Appendix. The whole procedure consists of four steps.

1. On each small stencil , , which contains and , there is a unique polynomial of degree at most which interpolates at the nodes in . Then we can get the integral

 JLi,r=α∫xixi−1e−α(xi−y)pr(y)dx=k∑j=0c(r)i−r−1+jvi−r−1+j, (2.18)

where the coefficients depend on and the cell size , but not on .

2. Similarly, on the big stencil , there is a unique polynomial of degree at most interpolating at the nodes in . Then we have

 JLi=α∫xixi−1e−α(xi−y)p(y)dx=2k−1∑j=0ci−k+jvi−k+j, (2.19)

and linear weights , such that

 JLi=k−1∑r=0drJLi,r, (2.20)

where .

3. Change the linear weights into nonlinear weights , using

 ωr=~ωr∑k−1s=0~ωs,  r=0,…,k−1, (2.21)

with

 ~ωr=dr(ϵ+βr)2.

Here, is a small number to avoid the denominator becoming zero, and we take in our numerical tests. The smoothness indicator , measuring the relative smoothness of the function in the stencil , is defined by

 βr=k∑l=1∫xixi−1Δx2l−1(∂lpr(x)∂xl)2dx, (2.22)
4. Lastly, we have

 JLi=k−1∑r=0ωrJLi,r. (2.23)

The process to obtain is mirror-symmetric to that for , with respect to the point .

When dealing with a Dirichlet or a Neumann boundary condition, we further need to extrapolate point values at the ghost points and using the point values in the domain, where

 x−i=x0−iΔx,  xM+i=xM+iΔx.

Following the idea in [40], a WENO extrapolation strategy is designed to get more robust results. For instance, when approximating the left ghost point , we consider the big stencil :

1. On each small stencil , , we can construct a polynomial of degree at most interpolating at the nodes in . And can be extrapolated by

 v−i=2k−2∑r=0drpr(x−i), (2.24)

where , , , , and .

2. Change the weights into nonlinear weights

 ωr=~ωr∑2k−2s=0~ωs,  r=0,…,2k−2, (2.25)

with

 ~ωr=dr(ϵ+βr)2.

Again, we take in our numerical tests to avoid the denominator becoming zero. And the smoothness indicator is defined by

 β0=Δx2, βr=r∑l=1∫x0x−1Δx2l−1(∂lpr(x)∂xl)2dx,  r=1,…,2k−2.
3. Finally, we have

 v−i=2k−2∑r=0ωrpr(x−i).

Similar, we can get the WENO extrapolation on right ghost points based on .

### 2.3 High order time discretization

Note that scheme (2.3) is first order accurate in time. Below, we will extend our schemes to high order accuracy via the implicit SSP RK method [21, 25].

To approximate the solutions of ordinary differential equation

 ut=F(u), (2.26)

an -stage and -th order RK method, denoted as , is defined as

 u(i)=un+Δts∑j=1aijF(u(j),tn+cjΔt),   i≤j≤s, (2.27a) un+1=un+Δts∑j=1bjF(u(j),tn+cjΔt), (2.27b)

with the standard assumption

 ci=s∑j=1aij.

Denote the matrix and the vector . In our scheme, we will use the diagonally implicit (DIRK) methods, for which A is lower triangular. Denote , , , and . Then, we can get . Hence

 Urk−Λ⋅F=un⋅Λ⋅A−1⋅1+(A−Λ)⋅A−1%Urk (2.28a) un+1=un+bT⋅F=un+bT⋅A−1(Urk−un⋅1) (2.28b)

Note that A is lower triangular. So is also lower triangular, with all diagnal elements equal to zero.

For problem (2.2), . Corresponding, each stage of (2.28a) can be rewritten as

 u(i)x+sgn(c⋅aii)αiu(i)=~αiun+i−1∑j=1~αiju(j), (2.29)

where . We can solve (2.29) with MOL method mentioned above. And finally, is obtained by (2.28b).

When problem (2.2) combining with the periodic boundary condition, the global coefficients or on each stage can be obtained by (2.16) or (2.17). For a Dirichlet boundary condition, the scheme may lose accuracy if we define or . To avoid order reduction, we follow the idea in [1]. In particular, denote for . To achieve -th order of accuracy, boundary values at intermediate stages are defined as

 Ark= k−2∑l=0(Al1⊗Δtl%cl)dldtlg1(tn)+(Ak−1⊗Δtk−1ck−1)dk−1dtk−1Gn+c1 (2.30)

for , where

 Ark=[Ark,1,…,Ark,s]T, Gn+c1=[g1(tn+c1Δt),…,g1(tn+csΔt)]T,

or

 Brk= k−2∑l=0(Al1⊗Δtl%cl)dldtlg2(tn)+(Ak−1⊗Δtk−1ck−1)dk−1dtk−1Gn+c2 (2.31)

for , where

 Brk=[Brk,1,…,Brk,s]T, Gn+c2=[g2(tn+c1Δt),…,g2(tn+csΔt)]T.

The similar idea is used when dealing with a Neumann boundary condition.

In our numerical tests, we will use and schemes to obtain third order accuracy and fourth order accuracy in time, respectively. Corresponding coefficients are listed as appendix.

### 2.4 Positivity-preserving limiter

The exact solution of (2.2) satisfies the positivity preserving principle, meaning that if the initial condition , then the exact solution for any time . In practice, it is highly desirable to have numerical schemes also satisfy such a property on the discrete level. It is even more crucial for the Vlasov simulation, since negative values may trigger some nonphysical precess in plasmas [33]. Recently, Zhang el at. developed a high order positivity-preserving limiters for finite volume WENO methods and finite element DG methods, which is able to maintain high order accuracy of the original scheme and easy to implement [46, 47, 48, 49]. Later, Xu et al. [45, 26, 42, 43, 44] developed a parameter maximum principle preserving flux limiter for finite difference WENO schemes by limiting the high order flux towards a first order monotone flux. However, both limiters are exclusively designed for MOL type methods and hence cannot be directly applied in the MOL setting. Below, we will formulate a new positivity-preserving (PP) limiter suited for the proposed WENO MOL scheme, and also show that such a limiter is able to preserve positive while maintaining the original high order accuracy.

Assume that we have for all at time . After obtaining the numerical solutions , we can rewrite our scheme as a conservative form

 un+1i=uni−Δt/Δx(^fi+1/2−^fi−1/2),

in which the artificial numerical fluxes are recursively defined by

 ^fi+1/2={0,i=−1,^fi−1/2−Δx/Δt(un+1i−uni),i=0,…,M.

If we are able to modify the numerical flux into , such that

 unewi=uni−Δt/Δx(~fi+1/2−~fi−1/2)≥0,  i=0,…,M, (2.32)

then the modified solution at time level will be non-negative automatically. Roughly speaking, the main idea of the limiter is to enforce the negative solution value back to zero by decreasing the outflow flux. For instance, when , we will reset the numerical flux if . Details of the limiter are shown below.

Firstly, we consider the case and a Dirichlet boundary condition. Numerical flux is modified into , starting from to ,

 ~fi+1/2=⎧⎨⎩~fi−1/2+Δx/Δtuni,if un+1,1i<0,^fi+1/2,otherwise, (2.33)

where

 un+1,1i=uni−Δt/Δx(^fi+1/2−~fi−1/2).

and , due to the Dirichlet boundary condition. For a Nuemann boundary condition, we can also use (2.33) to modify fluxes. Once the modified fluxes are obtained, the positivity-preserving solution is computed by (2.32).

For periodic boundary conditions, however, we need an additional step in order to preserve the positivity and periodicity of the solution. In particular, we first modify the numerical flux into by (2.33) i.e., starting from into , with

 ^^fi+1/2=⎧⎨⎩^fi−1/2+Δx/Δtuni,if un+1,1i<0,^fi+1/2,otherwise,

where

 un+1,1i=uni−Δt/Δx(^fi+1/2−^^fi−1/2).

For periodicity, we further set . Note that since the flux at is modified, the solution at may become negative. Hence, we need screen the solution one more time by the following procedure. Starting , we compute

 un+1,2i=uni−Δt/Δx(^^fi+1/2−~fi−1/2).

If we let

 ~fi+1/2=~fi−1/2+Δx/Δtuni.

Otherwise, let and for any with , let Finally, we are able to obtain the positivity-preserving solution via (2.32). Note that must be less than and hence .

When , wind blows from right to left, and we will reset the numerical flux if . Similar idea can be used for different boundary conditions.

Below, we show that the PP-limiter can preserve positivity while maintaining the previous high order. In particular, we have the following proposition.

###### Proposition 2.1.

Consider the WENO MOLscheme is used for solving (2.2) and assume that the numerical solution is -th order accurate in space and -th order accurate in time, i.e., the global error

 en+1i=|un+1i−u(xi,tn+1)|=O(Δxp+Δtq),∀i,

where is the exact solution. If the proposed PP-limiter (2.32) with (2.33) is applied, then

 ~en+1i=|unewi−u(xi,tn+1)|=O(Δxp−1+Δx−1Δtq),∀i

for the worst case, when .

Proof. Without loss of generality, we only consider (2.2) with and a Dirichlet boundary condition is imposed for brevity. Let be the first solution value that is less than zero starting from , i.e., and Then we have and such that Recalling that and and hence Also note that the modification of flux

 ~fi+1/2−^fi+1/2=Δx/Δtun+1i=O(Δxp+Δtq)

is a negative high order quantity.

Now, we consider the effect of the modification of flux to the next point value in terms of accuracy. First, we need to compute

 un+1,1i+1=uni+1−Δt/Δx(^fi+3/2−~fi+1/2).

Several cases need to be taken into account.

Case (1) If Then, by (2.33), and We have

 |unewi+1−un+1i+1|=|Δt/Δx(~fi+1/2−^fi+1/2)|=O(Δxp+Δtq).

Hence,

 ~en+1j+1=|unewi+1−u(xi+1,tn+1)|≤|unewi+1−un+1i+1|+en+1i+1=O(Δxp+Δtq).

Case (2) If Then, by (2.33), and . There are two cases we need to treat differently. If , then by the same argument we used to treat the case at , we have and hence . If , then we have

 |un+1,1i+1−un+1i+1|=|Δt/Δx(~fi+1/2−^fi+1/2)|=O(Δxp+Δtq).

Further, since and , combining the above estimate, we have

 |un+1i+1|=O(Δxp+Δtq).

Consequently,

 ~en+1j+1=|unewi+1−u(xi+1,tn+1)|=|u(xi+1,tn+1)|≤|un+1i+1|+en+1j+1=O(Δxp+Δtq).

We also have the following modification of the flux

 ~fi+3/2−^fi+3/2=Δx/Δtun+1i+~fi+1/2−^fi+1/2=Δx/Δt(un+1i+1+un+1i)=O(Δxp+Δtq). (2.34)

For both cases, we have shown that i.e., the PP-limiter preserves the order of accuracy. Moreover, for case (1) since the flux is not changed, we can move on to consider the solution value at with the same procedure as the case at . The situation of case (2) is more complicated when considering the subsequent effect of the modification of flux (2.34). Again, we need to consider two cases.

Case (2.1) If Similar to case (1) above, the flux at should remain the same, i.e., . By (2.34)

 |unewi+2−un+1i+2|=|Δt/Δx(~fi+3/2−^fi+3/2)|=O(Δxp+Δtq),

and consequently,

 ~en+1i+2=O(Δxp+Δtq).

Case (2.2) If . We can repeat the analysis of case (2) and obtain that

 |un+1i+2|=O(Δxp+Δtq),

which gives the estimate

 ~en+1i+2=O(Δxp+Δtq),

and the modification of the flux

 ~fi+5/2−^fi+5/2=Δx/Δt(un+1i+2+un+1i+1+un+1i)=O(Δxp+Δtq). (2.35)

More generally, assume that the fluxes are successively modified with at most constant many times, i.e., there exists a constant independent of and , such that if the fluxes for are successively changed based on

 ~fi+k+1/2=^fi+k+1/2+Δ