# Exponential DG methods for Vlasov equations

Nicolas Crouseilles Université de Rennes, Inria Rennes (Mingus team) and IRMAR
UMR CNRS 6625, F-35042 Rennes, France & ENS Rennes  Xue Hong Université de
Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042 Rennes,
France

###### Abstract

In this work, an exponential Discontinuous Galerkin (DG) method is proposed to
solve numerically Vlasov type equations. The DG method is used for space
discretization which is combined exponential Lawson Runge-Kutta method for
time discretization to get high order accuracy in time and space. In addition
to get high order accuracy in time, the use of Lawson methods enables to
overcome the stringent condition on the time step induced by the linear part
of the system. Moreover, it can be proved that a discrete Poisson equation is
preserved. Numerical results on Vlasov-Poisson and Vlasov Maxwell equations
are presented to illustrate the good behavior of the exponential DG method.

###### Contents

  1. 1 Introduction
  2. 2 Exponential DG method for 1D transport equation
  3. 3 Semi-discretization of some Vlasov models with Discontinuous Galerkin method
  4. 4 Vlasov-Maxwell 2dx-2dv
  5. 5 Numerical experiments
  6. 6 Conclusion
  7. A Appendix A: The DG matrix construction
  8. B Appendix B: Stability and error estimate for semi-discrete DG scheme.
  9. C Appendix C: Ker$(A)$ and projection $\Pi$
  10. D Appendix D: Lawson-Fourier method for Vlasov-Maxwell 1dx-2dv.

##  1 Introduction

In this work, we are interested in the numerical simulation of Vlasov type
equations using Eulerian based methods. Numerical approximation of Vlasov
equations has been the subject of a lot works since these models are widely
used to describe the dynamics of charged particles in a plasma through a
distribution function $f(t,x,v)$ with the time $t\geq 0$, $x$ the spatial
variables and $v$ the velocity variable. Hence, there have been a lot of
numerical methods that have been proposed to numerically solve Vlasov
equations. The so-called Particle-In-Cell(PIC) methods [30, 21] in which the
unknown is approximated by a sum of Dirac masses with a position and velocity
(macro-particles) that solves a differential system. Even if these methods are
efficient in high dimensions since only a spatial grid is required, they
however suffer from some numerical noise which make them hardly get an
accurate approximation. Indeed, the error slowly decreases when the number of
macro-particles increases, which turns out to be a drawback in low density
plasma region. On the other side, another family of methods have been
developed which uses a grid of the phase space $(x,v)$ like spectral methods
[20, 16] or finite differences/volumes methods [2, 1]. These methods enable to
get high order accurate approximation and as such, can capture fine physical
phenomena like Landau damping or filamentation in Vlasov equations.

However, due to the phase-space grid, these methods are quite costly both in
terms of memory and CPU point of view, in particular when high dimensions are
considered. Moreover, their stability is controlled by the so-called CFL
condition which imposes a constraint on the time step depending of the phase
space mesh refinement, which makes them very costly in practice. To overcome
this drawback, semi-Lagrangian have been developed [26, 28, 15] or arbitrary
Lagrangian-Eulerian methods [8, 31, 5, 18], which allow extra large time
stepping sizes with stability by tracking solutions along their
characteristics. For high dimensional problems, since (high order)
interpolation techniques are required which leads to huge computational cost.
Moreover, reaching high order accuracy in time are quite complicated. One
possibility is to use splitting methods which enable to deal with simple
subproblems that can even be solved exactly. However, the number of stages
required to get high order in time become prohibitive (see [12]).

For Vlasov type equations, the linear part induces the most stringent CFL
condition since the electromagnetic fields (which induce the nonlinear part)
are typically one order of magnitude smaller than the one of the linear
advection part. Based on this observation, exponential time integrator have
been proposed in which the linear part is solved exactly, and as such do not
suffer from the stability condition induced by the linear part, whereas the
nonlinear part is solved explicitely. These methods are very popular in a
number of applications ([17] and references therein) and enables to derive
easily high order methods in time since they are often based on a high order
Runge-Kutta method. Regarding the use of such time integrators for Vlasov
equations, we can quote [14, 12, 4], but these works are based on Fourier
techniques in space to approximate the linear part, and despite its simplicity
and its spectral accuracy, Fourier methods are quite limited in terms of
applications (cartesian domains) and suffer from Gibbs phenomenon when non
periodic boundary conditions are considered.

In the present work, we focus on exponential type method combined with
Discontinuous-Galerkin (DG) method in space to approximate Vlasov type
equations. The DG method is a class of finite element methods, in which the
approximation space contains completely discontinuous, piecewise polynomials
or other basis functions. High order accuracy can be obtained and complex
geometries with boundary conditions can be handled (see DG review article
[11]) which is important for physically relevant problems; moreover, one
element only communicates with its immediate neighbors which is very important
for parallelization capability [3] but also, thanks to this local data
structure, the matrices used in these methods are sparse which is an important
point if one wants to combine DG methods in space with exponential methods in
time.

Then, after a finite differences approximation of the Vlasov equation, a DG
method is employed for the space approximation using central fluxes. Indeed,
this choice is motivated by the fact that the DG matrix has a pure imaginary
spectrum which is not the case when monotone fluxes are considered. Moreover,
using central fluxes makes the treatment more easier compared to upwind flux.
Indeed, the latter requires to split the flux of the linear part into two
parts according to the sign of velocity $v$, which prevents a discrete Poisson
equation from being satisfied because of a lack consistency in the nonlinear
part. Once the semi-discrete equations is obtained, an ODE system has to be
solved in time. To do so, exponential time integrators are used to overcome
the stringent condition coming from the linear part, but the exponential of a
large matrix has to be computed. Thanks to a one side coupling between the
distribution function and the electromagnetic fields, it is possible to
compute explicitely the exponential of the matrix and to derive an efficient
numerical scheme which is high order in time, space and velocity, preserves
the total mass and a discrete Poisson equation. Some numerical illustrates the
good behavior of the method. In particular, the expected order are recovered
on a two-dimensional linear advection and a good agreement is obtained when we
compare the DG exponential method to the Fourier exponential method.

The rest of the paper is organized as follows. In Section 2, we present the
exponential DG method for one-dimensional (1D) linear transport problems.
Section 3 is dedicated to the construction of the exponential DG method for
Vlasov equations includes Vlasov-Ampère (1dx-1dv) and Vlasov Maxwell (1dx-2dv)
equations. In Section 4, we discuss some extensions of the exponential method
to high dimensional Vlasov-Maxwell equations. In Section 5, the capability of
the proposed exponential DG method is illustrated through several numerical
tests. Finally, after some concluding remarks, several appendices details some
specific aspects of the method.

##  2 Exponential DG method for 1D transport equation

We firstly consider the 1D transport equation:

$\left\\{\begin{aligned} &u_{t}+au_{x}=0,\quad x\in[x_{a},x_{b}],\\\
&u(0,x)=u_{0}(x).\end{aligned}\right.$ (1)

For simplicity, we assume periodic boundary conditions, and the velocity field
$a$ is a constant. Here we take $a$ as 1 for simplicity. We perform a
partition of the computational domain
$x_{a}=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots<x_{N+\frac{1}{2}}=x_{b}$ as the
mesh partition. Let $I_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ denote an
element of length $\Delta x_{j}=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$ and
define $\Delta x=\max_{j}\Delta x_{j}.$ For simplicity, we consider the
uniform mesh in this paper with $\Delta x_{j}=\Delta x=(x_{b}-x_{a})/N$. We
define the finite dimensional approximation space,
$V_{h}^{k}=\\{v_{h}:v_{h}|_{I_{j}}\in P^{k}(I_{j})\\}$, where $P^{k}(I_{j})$
denotes the set of polynomials of degree at most $k$ on $I_{j}$. For any
$\psi\in V_{h}^{k}$, we also denote the left limit of $\psi$ at cell boundary
as $\psi^{-}$ and the right limit as $\psi^{+}$. Multiply (1) by the test
function $\psi\in V_{h}^{k}$, integrate on cell $I_{j}$ and integrate by
parts, we end up with the semi-discrete DG scheme: find $u_{h}\in V_{h}^{k}$
such that

$\int_{{I}_{j}}(\partial_{t}u_{h}\psi)dx=-\hat{F}|_{{x}_{j+\frac{1}{2}}}\psi^{-}|_{{x}_{j+\frac{1}{2}}}+\hat{F}|_{{x}_{j-\frac{1}{2}}}\psi^{+}|_{{x}_{j-\frac{1}{2}}}+\int_{{I}_{j}}F\psi_{x}dx,\
j=1,...,N,$ (2)

where $F(u)\doteq u$ and $\hat{F}$ is chosen as either a central or upwind
fluxes

$\mbox{central flux:
}\;\hat{F}|_{{x}_{j\pm\frac{1}{2}}}=\frac{u^{-}+u^{+}}{2}|_{{x}_{j\pm\frac{1}{2}}},\;\;\mbox{upwind
flux: }\;\hat{F}|_{{x}_{j\pm\frac{1}{2}}}=u^{-}|_{{x}_{j\pm\frac{1}{2}}},$ (3)

Next, we consider $\xi_{j}^{m}\;(m=0,1,\dots,k)$ a basis of $P^{k}(I_{j})$ and
we choose a modal basis defined as $\xi^{m}_{j}(x)=((x-x_{j})/\Delta x)^{m}$
so that we have the representation
$u_{h}(t,x)|_{I_{j}}=\sum_{m=0}^{k}u_{j}^{m}(t)\xi_{j}^{m}(x)$, with
$u_{j}^{m}(t)$ the degree of freedom. The semi-discrete DG scheme can
eventually be written as an ordinary differential equation (ODE) satisfied by
the DG degrees of freedom $u^{m}_{j}(t)$ for $m=0,\dots,k$ and $j=1,\dots,N$.
Introducing the vector $\mathbf{u}(t)\in\mathbb{R}^{(k+1)N}$

$\mathbf{u}(t)=(u_{1}^{0},u_{1}^{1},\dots,u_{1}^{k},\;u_{2}^{0},u_{2}^{1},\dots,u_{2}^{k},\dots,u_{N}^{0},u_{N}^{1},\dots,u_{N}^{k})^{T}(t),$
(4)

the semi-discretized problem simply becomes

$\frac{d\mathbf{u}}{dt}=A\mathbf{u},\mbox{ with
}A\in\mathbb{M}_{{(k+1)N,(k+1)N}}(\mathbb{R}).$ (5)

The ’DG-matrix’ $A$ contains the DG approximation (2) for which the details
are given in Appendix A. From this semi-discrete in space formulation, a
Runge-Kutta discretization is classically used to get high order accuracy in
time [11, 33, 34]. But it seems also natural to use an exponential method for
time discretization, which turns out to be exact in this simple linear
transport case. Denoting $\mathbf{u}^{n}\approx\mathbf{u}(t^{n})$ the fully
discretized unknown (with $t^{n}=n\Delta t,n\in\mathbb{N}$ and $\Delta t>0$
the time step), the exponential-DG scheme thus writes

$\mathbf{u}^{n+1}=\exp(A\Delta t)\mathbf{u}^{n},\;\;\;\forall n\in\mathbb{N},$
(6)

with $\mathbf{u}^{0}=\mathbf{u}_{0}$ ($\mathbf{u}_{0}$ being the degrees of
freedom of the initial condition $u_{0}$ in (1)). Note that from
$\mathbf{u}^{n}\in\mathbb{R}^{(k+1)N}$ (whose components are denoted by
$(u_{j}^{m})^{n},j=1,\dots,N$ and $m=0,\dots,k$, following (4)), it is
possible to reconstruct a piecewise polynomial function $u^{n}_{h}\in
V_{h}^{k}$ from

$u_{h}^{n}(x)=\sum_{m=0}^{k}(u_{j}^{m})^{n}\xi_{j}^{m}(x),\;\;\forall x\in
I_{j},\;\;j=1,\dots,N.$ (7)

The properties of such an approximation obviously depends on the structure of
the DG-matrix $A$ which is discussed now. From the calculations (given in
Appendix A), the DG-matrix $A$ obtained with the central flux and periodic
boundary conditions enjoys a circulant tri-diagonal block structure so that it
can be written as:

$A=\frac{1}{\Delta
x}\left(\begin{array}[]{lllll}C_{1}&C_{2}&\bf{0}&\ldots&C_{3}\\\
C_{3}&C_{1}&C_{2}&\bf{0}&\ldots\\\ \bf{0}&\ddots&\ddots&\ddots&C_{2}\\\
C_{2}&\bf{0}&\ldots&C_{3}&C_{1}\end{array}\right),\;\;\mbox{ with
}C_{j}=M^{-1}D_{j}\;(j=1,2,3)$ (8)

where the matrix elements of the matrices
$M,D_{j}\in\mathbb{M}_{{(k+1),(k+1)}}(\mathbb{R})$ are given by (for
$\ell,m=1,\dots,k+1$, see Appendix A for details)

$\displaystyle
M_{\ell,m}=\frac{(1/2)^{m+\ell-1}}{m+\ell-1}[1-(-1)^{m+\ell-1}],\\!\\!\\!\\!\\!\\!$
$\displaystyle(D_{2})_{\ell,m}=(-1)^{m}(1/2)^{m+\ell-1},$
$\displaystyle(D_{1})_{\ell,m}=(1/2)^{m+\ell-2}\Big{(}\frac{\ell-1}{m+\ell-2}-\frac{1}{2}\Big{)}[1-(-1)^{m+\ell-2}],\\!\\!\\!\\!\\!\\!$
$\displaystyle(D_{3})_{\ell,m}=(-1)^{\ell-1}(1/2)^{m+\ell-1},$ (9)

with $(D_{1})_{1,1}=0$ by convention. The choice of central fluxes implies the
matrix $A$ is diagonalizable and the eigenvalues are pure imaginary. This has
been checked numerically and some discussions are performed in the following
remarks.

######  Remark 1.

In [29], the author proposes a way to deduce the eigenvalues of
$A\in\mathbb{M}_{{(k+1)N,(k+1)N}}(\mathbb{R})$ from the eigenvalues of some
matrices of size $(k+1)$, which can be computed explicitly for small values of
$k$ (numerically for larger $k$). Considering $(\rho_{j})_{j=0,\dots,N-1}$ the
$N$-th roots of the unity ($\rho_{j}^{N}=1$ for $j=0,\dots,N-1$), the $(k+1)N$
eigenvalues of $A$ given by (8) can be deduced from the $(k+1)$ eigenvalues of
${\cal C}_{j}=C_{1}+\rho_{j}C_{2}+\rho_{j}^{N-1}C_{3}$ for $j=0,\dots,N-1$,
where ${\cal C}_{j}\in\mathbb{M}_{{(k+1),(k+1)}}$. Then, we checked
numerically that the eigenvalues of ${\cal C}_{j}$ are pure imaginary for all
$j=0,\dots,N-1$, and we deduce from [29] that it is also true for the
eigenvalues of $A$.

######  Remark 2.

We explore another way to check the eigenvalues of $A$ given by (8) are pure
imaginary by using symbolic software. Denoting $P_{A}(\lambda)$ the
characteristic polynomial of $A$, we made the following observations

  * •

odd case: $(k+1)N=2d+1$. In this case, we have
$P_{A}(\lambda)=\lambda\sum_{\ell=0}^{d}a_{2\ell}\lambda^{2\ell}$ with
$a_{2\ell}\in\mathbb{R}$ and the roots can be written as
$\left\\{0,\lambda_{j},\bar{\lambda}_{j},j=1,\dots,d\right\\}$, in particular
$0$ is a simple eigenvalue in this case. Since
$P_{A}(-\lambda)=P_{A}(\lambda)$, we deduce Re$(\lambda_{j})=0$.

  * •

even case: $(k+1)N=2d$. In this case, we have
$P(\lambda)=\lambda^{2}\sum_{\ell=0}^{d-1}a_{2\ell}\lambda^{2\ell}$ with
$a_{2\ell}\in\mathbb{R}$ and the roots can be written as
$\left\\{0,\lambda_{j},\bar{\lambda}_{j},j=1,\dots,d-1\right\\}$, in
particular $0$ is a double eigenvalue in this case. Since
$P_{A}(\lambda)=P_{A}(-\lambda)$, we deduce Re$(\lambda_{j})=0$.

We can now study the stability of the numerical scheme. To do so, we write the
following proposition.

######  Proposition 1.

Let consider the matrix $A\in\mathbb{M}_{(k+1)N,(k+1)N}(\mathbb{R})$ given by
(8)-(9) There exists $C>0$ such that, for any time $t$ and any $k,N$, we have
$\|\exp(At)\|\leq C$, with $\|\cdot\|$ an induced matrix norm.

###### Proof.

First, we write $A=\Delta x^{-1}A_{1}$ (with $\Delta x=(x_{b}-x_{a})/N$) and
since $A$ is diagonalizable, there exist $P$ invertible and $D$ diagonal such
that $A=\Delta x^{-1}PDP^{-1}$. Let us remark from (8)-(9) that $A_{1}$ does
not depend on the space mesh $\Delta x$ (and then dos not depend on $N$), so
does the matrix $P$. Thus, there exists $C>0$ (independent of $N$) such that
cond$(P)\equiv\|P\|\|P^{-1}\|\leq C$, where $\|\cdot\|$ denotes an induced
matrix norm. Now, since the eigenvalues of $A$ are pure imaginary for all
$k,N$, we have
$D_{j,j}=i\lambda_{j},\lambda_{j}\in\mathbb{R},j=1,2,\dots,(k+1)N$. Finally,
we get for all $t$

$\|\exp(At)\|=\|P\exp(\Delta x^{-1}D)P^{-1}\|\leq\mbox{cond}(P)\;\|\exp(\Delta
x^{-1}D)\|=\mbox{cond}(P)\leq C.$

∎

We end this section by proving an error estimate for exponential DG method.

######  Proposition 2.

Let $u(t,x)$ the exact solution of 1D transport problem (1) with a smooth
initial condition $u_{0}$ and let $u_{h}^{n}\in V_{h}^{k}$ the numerical
solution, $n=0,\dots,N$ (with $N=T/\Delta t$, $T$ being the final time and
$\Delta t$ the time step) obtained from (6)-(7) where $A$ is the DG-matrix
given by (8)-(9). Then we have the following $L^{2}$-norm error estimate:

$||u(t^{n},\cdot)-u_{h}^{n}||_{L^{2}}\leq C\Delta x^{k},\;\;\;$

###### Proof.

First, we introduce $u_{h}(t,x)$ the exact solution of the semi-discrete DG
scheme (5). The classical DG projection analysis gives
$||u(t^{n},\cdot)-u_{h}(t^{n},\cdot)||_{L^{2}}\leq C\Delta x^{k},$ where $C$
is a positive constant independent on $\Delta x$ (see the details in [27, 22,
19, 25, 33, 34], but a simple proof of stability and error estimate is given
in Appendix B). Since the exponential method exactly solves the semi-discrete
DG scheme, we have $||u_{h}(t^{n},\cdot)-u_{h}^{n}||_{L^{2}}=0$. Finally we
have

$||u(t^{n},\cdot)-u_{h}^{n}||_{L^{2}}\leq||u(t^{n},\cdot)-u_{h}(t^{n},\cdot)||_{L^{2}}+||u_{h}(t^{n},\cdot)-u_{h}^{n}||_{L^{2}}\leq
C\Delta x^{k}.$

∎

######  Remark 3.

It is worthy to be mentioned that the $(k+1)$th order optimal convergence rate
have been proved for DG with monotone flux (see [19, 25, 33, 34]). However,
$k$th order sub-optimal convergence rate is proved for DG with central flux in
[23] and a discussion is performed according to the oddness of $k$.

##  3 Semi-discretization of some Vlasov models with Discontinuous Galerkin
method

In this section, we consider the numerical approximation of Vlasov-Maxwell
equations using DG framework in space (as presented in the previous section)
and finite differences in the velocity direction. The semi-discretization (in
both space and velocity) is presented and we will see that the so-obtained ODE
system is amenable to Lawson time integrators.

We first present the methodology on the $1d_{x}-1d_{v}$ case on the Vlasov-
Ampère system and then we consider the $1d_{x}-2d_{v}$ case on the Vlasov-
Maxwell system.

###  3.1 Vlasov-Ampère equation

The equation we address is the following Vlasov-Ampère model satisfied by the
distribution function $f(t,{x},{v})\geq 0$ and the electric field
$E(t,x)\in\mathbb{R}$ with $t\geq 0,x\in[0,L]\;(L>0)$ and $v\in\mathbb{R}$,

$\left\\{\begin{aligned} &\frac{\partial f}{\partial t}+{v}\frac{\partial
f}{\partial{x}}+E\frac{\partial f}{\partial{v}}=0,\\\
&\partial_{t}E=-\int_{\mathbb{R}}vf\mathrm{d}v+\bar{J},\;\;\;\;\bar{J}=\frac{1}{L}\int_{0}^{L}\int_{\mathbb{R}}vfdxdv,\\\
\end{aligned}\right.$ (10)

with the initial conditions $(f_{0}(x,v),E_{0}(x))$ such that the Poisson
equation is satisfied initially
$\partial_{x}E_{0}=\int_{\mathbb{R}}f_{0}dv-\bar{\rho}$, with
$\bar{\rho}=(1/L)\int_{0}^{L}\int_{\mathbb{R}}f_{0}dvdx$ and periodic boundary
conditions are imposed in space. The Vlasov-Ampère system is equivalent the
Vlasov-Poisson model where the electric field satisfies the Poisson equation
$\partial_{x}E=\int_{\mathbb{R}}f{d}{v}-\bar{\rho}$.

####  3.1.1 Semi-discretization

We shall use a DG method in the space direction $x$ as presented in the
previous section and we consider a truncated domain $[-v_{\max},v_{\max}]$ in
the velocity direction discretized by $v_{j}=-v_{\max}+j\Delta
v,j=0,\dots,N_{v}$, $\Delta v=2v_{\max}/N_{v}$ being the velocity mesh step.
We firstly present DG discretization for ${E}$ in the $x$-direction with
$I_{i}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}],\ i=1,...,N_{x}$ ($N_{x}$ being
the number of cells):

$E(t,x)\approx
E_{h}(t,x)=\sum_{i=1}^{N_{x}}E_{h}(t,x)|_{I_{i}}=\sum_{i=1}^{N_{x}}\sum_{m=0}^{k}{E}_{i}^{m}(t)\xi^{m}(x).$

Then, considering finite difference method for $f$ in $v$ direction, we
consider the DG approximation in $x$-direction through

$f(t,x,v_{j})\approx
f_{h}(t,x,v_{j})=\sum_{i=1}^{N_{x}}\sum_{m=0}^{k}{f}_{i}^{m}(t,v_{j})\xi^{m}(x),\
j=0,...,N_{v}.$

As in the previous section, we denote
$\mathbf{f}_{j}(t)\in\mathbb{R}^{(k+1)N_{x}}$ the vector of the DG
coefficients $f^{m}_{i}(t,v_{j})$ of $f_{h}(t,x,v_{j})$ using DG in space and
evaluated at the velocity grid $v_{j}$ whereas
$\mathbf{E}\in\mathbb{R}^{(k+1)N_{x}}$ denotes the vector of DG coefficients
$E^{m}_{i}$ of $E_{h}(t,x)$.

For Vlasov-Ampère equation (10), we have the following DG scheme with the DG
representation of $E$ and $f$:

$\displaystyle\sum_{m=0}^{k}\left[\Big{(}\partial_{t}{f}_{i}^{m}(t,v_{j})\xi^{m},\xi^{\ell}\Big{)}_{I_{i}}-v_{j}\Big{(}{f}_{i}^{m}(t,v_{j})\xi^{m},\partial_{x}\xi^{\ell}(x)\Big{)}_{I_{i}}\right]+v_{j}\left[\\{f_{h}(t,x,v_{j})\\}\xi^{\ell}\right]^{i+1/2}_{i-1/2}$
(11)
$\displaystyle+\sum_{m=0}^{k}\Big{(}\sum_{n=0}^{k}E_{i}^{n}\xi^{n}({\mathcal{D}f_{i}^{m}})(t,v_{j})\xi^{m},\xi^{\ell}\Big{)}_{I_{i}}=0,$

where we used the central flux $\\{f_{h}(t,x,v_{j})\\}|_{x_{i\pm
1/2}}=\frac{1}{2}(f_{h}(t,x_{i\pm 1/2}^{+},v_{j})+f_{h}(t,x_{i\pm
1/2}^{-},v_{j}))$, $\ell=0,1,2,...,k,\ i=1,2,...N_{x}$ and
${\mathcal{D}f}(v_{j})$ denotes a discrete approximation of
$(\partial_{v}f)(v_{j})$ (an example would be
$({\mathcal{D}f})(v_{j})=\frac{f(v_{j+1})-f(v_{j-1})}{2\Delta v}$ but any
higher order finite difference approximation can be used). We denote

$\mathbf{f}_{j,i}(t)=({f}_{i}^{0}(t,v_{j}),{f}_{i}^{1}(t,v_{j}),....,{f}_{i}^{k}(t,v_{j}))^{T},\
i=1,2,...,N_{x},$
$\mathbf{E}_{i}(t)=({E}_{i}^{0}(t),{E}_{i}^{1}(t),....,{E}_{i}^{k}(t))^{T},\
i=1,2,...,N_{x},$

and

$\mathbf{{\mathcal{D}f}}_{j,i}(t)=(({\mathcal{D}f}_{i}^{0})(t,v_{j}),({\mathcal{D}f}_{i}^{1})(t,v_{j}),....,({\mathcal{D}f}_{i}^{k})(t,v_{j}))^{T},\
i=1,2,...,N_{x}.$

We can rewrite the DG discretization as an ODE system of size $(k+1)N_{x}$ for
each $j=1,\dots,N_{v}$

$\displaystyle\Delta x\left(\begin{array}[]{llll}M&&&\\\ &M&&\\\ &&\ddots&\\\
&&&M\\\
\end{array}\right)\frac{d}{dt}\left(\begin{array}[]{llll}\mathbf{f}_{j,1}\\\
\mathbf{f}_{j,2}\\\ \vdots\\\ \mathbf{f}_{j,N_{x}}\end{array}\right)$
$\displaystyle-v_{j}\left(\begin{array}[]{llll}D_{1}&D_{2}&\ldots&D_{3}\\\
D_{3}&D_{1}&D_{2}&\\\ &\ddots&\ddots&D_{2}\\\ D_{2}&\ldots&D_{3}&D_{1}\\\
\end{array}\right)\left(\begin{array}[]{llll}\mathbf{f}_{j,1}\\\
\mathbf{f}_{j,2}\\\ \vdots\\\ \mathbf{f}_{j,N_{x}}\end{array}\right)$ (12)
$\displaystyle+\left(\begin{array}[]{llll}B_{1}&&&\\\ &B_{2}&&\\\ &&\ddots&\\\
&&&B_{N_{x}}\\\
\end{array}\right)\left(\begin{array}[]{llll}\mathbf{({\mathcal{D}f})}_{j,1}\\\
\mathbf{({\mathcal{D}f})}_{j,2}\\\ \vdots\\\
\mathbf{({\mathcal{D}f})}_{j,N_{x}}\end{array}\right)={\bf 0},$

where the matrices $M,D_{1},D_{2},D_{3}$ are the same as in the previous
section (see also Appendix A) and $B_{i}(i=1,\dots,N_{x})$ are matrices of
size $k+1$ with elements
$(B_{i})_{\ell,m}=(\sum_{n=0}^{k}E_{i}^{n}\xi^{n}\xi^{m},\xi^{\ell})_{I_{i}}$.

Introducing now the following vector containing the degrees of freedom of
$f_{h}$ and $E_{h}$

$\displaystyle\mathbf{f}_{j}(t)\\!\\!$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\\!\\!({f}_{1}^{0}(t,v_{j}),\dots,{f}_{1}^{k}(t,v_{j}),{f}_{2}^{0}(t,v_{j}),\dots,{f}_{2}^{k}(t,v_{j}),\dots,{f}_{N_{x}}^{0}(t,v_{j}),\dots,{f}_{N_{x}}^{k}(t,v_{j}))^{T},$
$\displaystyle\mathbf{E}(t)\\!\\!$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\\!\\!({E}_{1}^{0}(t),\dots,{E}_{1}^{k}(t),{E}_{2}^{0}(t),\dots,{E}_{2}^{k}(t),\dots,{E}_{N_{x}}^{0}(t),\dots,{E}_{N_{x}}^{k}(t))^{T},$

and $\mathbf{{\mathcal{D}f}}_{j}$ is defined similarly as $\mathbf{f}_{j}$, we
can rewrite the DG scheme (12) as

$\partial_{t}\mathbf{f}_{j}=v_{j}A\mathbf{f}_{j}-\tilde{\mathbf{E}}(\mathbf{{\mathcal{D}f}})_{j},$
(13)

$A\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}$ is the DG-matrix (8) and
$\tilde{\mathbf{E}}\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}$ is a block diagonal
matrix composed of $N_{x}$ block matrices of size $(k+1)\times(k+1)$ defined
by $(\Delta xM)^{-1}B_{i},i=1,\dots,N_{x}$. Let remark that we consider
central finite differences method to approximate
$(\mathbf{{\mathcal{D}f}})_{j}$ to avoid to discuss the sign of matrix
$\tilde{\mathbf{E}}\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}$ compared with upwind
FD method.

Let us now discuss the discretization of the Ampère equation. Since our goal
is to find a consistent discretization that is compatible with a discrete
Poisson equation, we first discuss how to solve the Poisson equation. Using
the above discretization, we will use the DG matrix $A$ which is an
approximation of $(-\partial_{x})$. Thus, a direct approximation of the
initial Poisson equation
$\partial_{x}E(0,x)=\int_{\mathbb{R}}f_{0}(x,v){d}{v}-\bar{\rho}$ would be
$-A\mathbf{E}^{0}=\sum_{j}\mathbf{f}^{0}_{j}\Delta v-\bar{\rho}$ (with
$\mathbf{E}^{0}$ and $\mathbf{f}^{0}_{j}$ the degrees of freedom of $E(0,x)$
and $f_{0}(x,v_{j})$. However, as mentioned in Remark 2,
$A\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}$ is not invertible and we then
introduce $\Pi\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}$ the projection onto the
Ker$(A)$ so that $(A+\Pi)$ is invertible on $R(A)$ with $R(A)$ the range of
$A$. Here we impose condition $\Pi\mathbf{E}^{0}=0$ to preserve the uniqueness
of the solution $\mathbf{E}\in\mathbb{R}^{(k+1)N_{x}}$, which is similarly as
the constraint $\int_{\mathbb{R}}E(x,v){d}{x}=0$ for Poisson equation itself.
We then consider the following discretized Poisson equation

$-(A+\Pi)\mathbf{E}^{0}=({\bf{1}}-\Pi)\Big{(}\sum_{j}\mathbf{f}^{0}_{j}\Delta
v-\bar{\rho}\Big{)}=({\bf{1}}-\Pi)\sum_{j}\mathbf{f}^{0}_{j}\Delta v,$ (14)

with ${\bf{1}}$ the identity matrix of size $(k+1)N_{x}$ and where in the last
equality, we used the fact that constants belong to Ker$(A)$ (see Appendix C
for details).

We deduce the discretization of Ampère equation from the time derivative of
the discretized Poisson equation inspired from (14), that is:
$-(A+\Pi)\mathbf{E}(t)=({\bf{1}}-\Pi)\sum_{j}\mathbf{f}_{j}(t)\Delta v$.
Considering the time derivative of the latter equation and using (13) leads to

$\displaystyle-(A+\Pi)\partial_{t}\mathbf{E}(t)$ $\displaystyle=$
$\displaystyle({\bf{1}}-\Pi)\partial_{t}\sum_{j}\mathbf{f}_{j}(t)\Delta
v=({\bf{1}}-\Pi)\sum_{j}v_{j}A\mathbf{f}_{j}(t)\Delta v$ $\displaystyle=$
$\displaystyle A({\bf{1}}-\Pi)\sum_{j}v_{j}\mathbf{f}_{j}(t)\Delta
v=(A+\Pi)({\bf{1}}-\Pi)\sum_{j}v_{j}\mathbf{f}_{j}(t)\Delta v,$

where we used $\sum_{j}{\cal D}{\bf f}_{j}=0$ and some relations between $A$
and $\Pi$. Hence, we consider the following DG discretization of the Ampère
equation

$\partial_{t}\mathbf{E}(t)=-({\bf{1}}-\Pi)\sum_{j}v_{j}\mathbf{f}_{j}(t)\Delta
v.$ (15)

Finally, gathering (13) and (15) enables to get the following semi-discretized
scheme for the Vlasov-Ampère system

$\left\\{\begin{aligned}
&\partial_{t}\mathbf{f}_{j}=v_{j}A\mathbf{f}_{j}-\tilde{\mathbf{E}}(\mathbf{{\mathcal{D}f}})_{j},\\\
&\partial_{t}\mathbf{E}=-\sum_{j}v_{j}\mathbf{f}_{j}\Delta
v+\sum_{j}v_{j}\Pi\mathbf{f}_{j}\Delta v.\end{aligned}\right.$ (16)

In view of the time discretization, we introduce the following vector of semi-
discrete unknown
$U=(\vec{\mathbf{f}},\mathbf{E})\in\mathbb{R}^{(k+1)N_{x}(N_{v}+1)}$ with
$\vec{\mathbf{f}}=(\mathbf{f}_{1},\mathbf{f}_{2},\dots,\mathbf{f}_{N_{v}})\in\mathbb{R}^{(k+1)N_{x}N_{v}}$.
Then the previous system (16) can be rewritten as

$\partial_{t}U=LU+N(U),$ (17)

with $L\in\mathbb{M}_{(k+1)N_{x}(N_{v}+1),(k+1)N_{x}(N_{v}+1)}$ given by

$\displaystyle L$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{lllllllccc}\;\;\;\;\;\;{v_{1}}A&\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}\;\;\;\;\ldots&\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&{\bf 0}_{\tilde{N},\tilde{N}}\\\ \;\;\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&\;\;\;\;{v_{2}}A&{\bf
0}_{\tilde{N},\tilde{N}}\;\;\;\;\ldots&\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&\vdots\\\
\;\;\;\;\;\;\vdots&\;\;\;\;\ddots&\ddots&\;\;\;\;\vdots&\vdots\\\
\;\;\;\;\;\;{\bf 0}_{\tilde{N},\tilde{N}}&\;\;\;\;\ldots&{\bf
0}_{\tilde{N},\tilde{N}}&\;\;\;\;{v_{N_{v}}}A&{\bf 0}_{\tilde{N},\tilde{N}}\\\
\\!\\!\\!-\Delta vv_{1}({\bf 1}-\Pi)&-\Delta vv_{2}({\bf
1}-\Pi)&\ldots&-\Delta vv_{N_{v}}({\bf 1}-\Pi)&{\bf
0}_{\tilde{N},\tilde{N}}\end{array}\right)$ (23)

where we denote $\tilde{N}=(k+1)N_{x}$ and ${\bf 1}={\bf
1}_{\tilde{N},\tilde{N}}$ the identity matrix of size
$\tilde{N}\times\tilde{N}$. Finally,
$U,N(U)\in\mathbb{R}^{(k+1)N_{x}(N_{v}+1)}$ are given by

$\displaystyle U$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{llll}{\mathbf{f}}_{1}\\\
{\mathbf{f}}_{2}\\\ \vdots\\\ \\!\\!\\!{\mathbf{f}}_{N_{v}}\\!\\!\\!\\!\\\
\mathbf{E}\end{array}\right),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;N(U)=\left(\begin{array}[]{llllccc}-\tilde{\mathbf{E}}(\mathbf{{\mathcal{D}f}})_{1}\\\
-\tilde{\mathbf{E}}(\mathbf{{\mathcal{D}f}})_{2}\\\ \;\;\;\;\vdots\\\
-\tilde{\mathbf{E}}(\mathbf{{\mathcal{D}f}})_{N_{v}}\\\ \;\;\;\;{\bf
0}\end{array}\right).$

####  3.1.2 Time discretisation

The goal of this part is to present time discretization of (17) to get a fully
discretized scheme of the Vlasov-Ampère system (10). The form (17) is amenable
to exponential scheme [17, 14] which is motivated by the fact that the linear
part acts on a different scale compared to the nonlinear part in Vlasov type
problems. Moreover, as discussed in Section 2, the linear part can be computed
exactly thanks to the exponential. Among the exponential schemes, we shall use
the Lawson class of methods for stability reasons [14].

Denoting $U^{n}=({\vec{\bf f}^{n},{\bf E}^{n}})\approx({\vec{\bf
f}(t^{n}),{\bf E}(t^{n})})=U(t^{n})$ with $t^{n}=n\Delta t,\forall
n\in\mathbb{N}$ ($\Delta t>0$ being the time step), the simplest (first order
in time) Lawson method can be written as

$U^{n+1}=\exp(\Delta tL)U^{n}+\Delta t\exp(\Delta tL)N(U^{n}).$ (25)

High order methods can be obtained from Runge-Kutta methods using the
corresponding Butcher tableau [14].

The key point of exponential methods lies in the computation of $\exp(\Delta
tL)$ with $L$ given by (23). To compute $\exp(At)\ \forall t>0$, one considers
the linear part only $\partial_{t}U=LU$. First, we observe that the
distribution function $\mathbf{f}_{j}$ part is decoupled from the electric
field part, so that it can be solved directly and we have, similarly as in
Section 2

$\mathbf{f}_{j}(t)=\exp(v_{j}At)\mathbf{f}_{j}(0).$ (26)

Now, let consider the electric field equation

${\partial_{t}\mathbf{E}}=-\sum_{j}v_{j}({\bf 1}-\Pi)\mathbf{f}_{j}(t)\Delta
v,$

which gives, after integrating it in time

$\mathbf{E}(t)=\mathbf{E}(0)-\sum_{j}\int_{0}^{t}v_{j}({\bf
1}-\Pi)\mathbf{f}_{j}(s)\Delta v\mathrm{d}s.$

Now, replacing $\mathbf{f}_{j}(s)$ by $\exp(v_{j}As)\mathbf{f}_{j}(0)$ from
(26) enables to get an explicit expression of $\mathbf{E}(t)$. Indeed, using
properties used to derive (14) we have

$\displaystyle\mathbf{E}(t)$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}\int_{0}^{t}v_{j}({\bf
1}-\Pi)\mathbf{f}_{j}(s)\Delta
v\mathrm{d}s=\mathbf{E}(0)-\sum_{j}\int_{0}^{t}v_{j}({\bf
1}-\Pi)\exp(v_{j}As)\mathbf{f}_{j}(0)\mathrm{d}s\Delta v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}\int_{0}^{t}v_{j}({\bf
1}-\Pi)\Bigg{[}\sum_{k=0}^{\infty}\frac{(v_{j}As)^{k}}{k!}\Big{(}\Pi\mathbf{f}_{j}(0)+({\bf
1}-\Pi)\mathbf{f}_{j}(0)\Big{)}\Bigg{]}\mathrm{d}s\Delta v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}v_{j}({\bf
1}-\Pi)\int_{0}^{t}\Bigg{[}\Pi\mathbf{f}_{j}(0)+\sum_{k=0}^{\infty}\frac{(v_{j}(A+\Pi)s)^{k}}{k!}({\bf
1}-\Pi)\mathbf{f}_{j}(0)\Bigg{]}\mathrm{d}s\Delta v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}v_{j}({\bf
1}-\Pi)\Bigg{[}t\Pi\mathbf{f}_{j}(0)+\sum_{k=0}^{\infty}\frac{(v_{j}(A+\Pi))^{k}t^{k+1}}{(k+1)!}({\bf
1}-\Pi)\mathbf{f}_{j}(0)\Bigg{]}\Delta v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}\Bigg{[}({\bf
1}-\Pi)\sum_{k=0}^{\infty}(A+\Pi)^{-1}\frac{(v_{j}(A+\Pi)t)^{k+1}}{(k+1)!}({\bf
1}-\Pi)\mathbf{f}_{j}(0)\Bigg{]}\Delta v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}({\bf
1}-\Pi)(A+\Pi)^{-1}(\exp(v_{j}At)-{\bf 1})({\bf 1}-\Pi)\mathbf{f}_{j}(0)\Delta
v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)-\sum_{j}(A+\Pi)^{-1}(\exp(v_{j}At)-{\bf 1})({\bf
1}-\Pi)\mathbf{f}_{j}(0)\Delta v$ $\displaystyle=$
$\displaystyle\mathbf{E}(0)+\Delta v\sum_{j}(A+\Pi)^{-1}({\bf
1}-\exp(v_{j}At))\mathbf{f}_{j}(0).$

Hence, denoting $A+\Pi$ as $\tilde{A}$, we deduce from the above calculation
the explicit expression of $\exp(Lt)$

$\exp(Lt)=\left(\begin{array}[]{lllllllccc}\;\;\;\;\;\;\;e^{v_{1}At}&\;\;\;\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}\;\;\;\;\;\;\;\ldots&\;\;\;\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&{\bf 0}_{\tilde{N},\tilde{N}}\\\ \;\;\;\;\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&\;\;\;\;\;\;\;e^{v_{2}At}&{\bf
0}_{\tilde{N},\tilde{N}}\;\;\;\;\;\;\;\ldots&\;\;\;\;\;\;\;{\bf
0}_{\tilde{N},\tilde{N}}&\vdots\\\
\;\;\;\;\;\;\;\;\vdots&\;\;\;\;\;\;\;\ddots&\ddots&\;\;\;\;\;\;\;\vdots&\vdots\\\
\;\;\;\;\;\;\;\;{\bf 0}_{\tilde{N},\tilde{N}}&\;\;\;\;\;\;\;\ldots&{\bf
0}_{\tilde{N},\tilde{N}}&\;\;\;\;\;\;\;e^{v_{N_{v}}At}&{\bf
0}_{\tilde{N},\tilde{N}}\\\ \\!\\!\Delta v\tilde{A}^{-1}({\bf
1}-e^{v_{1}At})&\Delta v\tilde{A}^{-1}({\bf 1}-e^{v_{2}A\Delta
t})&\ldots&\Delta v\tilde{A}^{-1}({\bf 1}-e^{v_{N_{v}}At})&{\bf
1}\end{array}\right).$ (27)

Hence, the exponential-DG scheme for the Vlasov-Ampère equation corresponds to
(25)-(27). For this scheme, one can prove in the following proposition that a
discrete Poisson equation is satisfied for each iteration.

######  Proposition 3.

The exponential DG method (25)-(27) (and its generalization to high order
Lawson Runge-Kutta) satisfied by $U^{n}=(\vec{\bf{f}}^{n},{\bf E}^{n})$
preserves the following discretized Poisson equation:

$(A+\Pi)\mathbf{E}^{n}=-\sum_{j}({\bf 1}-\Pi)\mathbf{f}^{n}_{j}\Delta
v,\;\;\;\forall n\in\mathbb{N}^{\star},$

provided that it is satisfied at the initial time $n=0$. Here, $A$ is the DG
matrix given by (8)-(9), $\Pi$ is the orthogonal projection onto Ker$(A)$ and
${\bf{1}}$ is the identity matrix of size $(k+1)N_{x}$,

###### Proof.

We present the proof for first order Lawson case (forward Euler), the proof
can be generalized to arbitrary explicit Runge-Kutta scheme. First, we assume
the Poisson equation $(A+\Pi)\mathbf{E}^{0}=-\sum_{j}({\bf
1}-\Pi)\mathbf{f}_{j}^{0}\Delta v$ holds at the initial time.
Next, from the scheme (25) with $\exp(Lt)$ given by (27), we have

$\mathbf{f}_{j}^{n+1}=\exp(v_{j}A\Delta t)\mathbf{f}_{j}^{n}-\Delta
t\exp(v_{j}A\Delta t)(\tilde{\mathbf{E}}\mathbf{{\mathcal{D}f}})_{j}^{n},$
(28)

whereas for the ${\bf E}$ component, we have

$\displaystyle\mathbf{E}^{n+1}$ $\displaystyle=$
$\displaystyle\mathbf{E}^{n}+\Delta v\sum_{j}(A+\Pi)^{-1}({\bf
1}-\exp(v_{j}A\Delta t))\mathbf{f}_{j}^{n}$ $\displaystyle-\Delta t\Delta
v\sum_{j}(A+\Pi)^{-1}({\bf 1}-\exp(v_{j}A\Delta
t))(\tilde{\mathbf{E}}\mathbf{{\mathcal{D}f}})_{j}^{n}$ $\displaystyle=$
$\displaystyle\mathbf{E}^{n}+\Delta v\sum_{j}(A+\Pi)^{-1}({\bf 1}-\Pi)({\bf
1}-\exp(v_{j}A\Delta t))\mathbf{f}_{j}^{n}$ $\displaystyle-\Delta t\Delta
v\sum_{j}(A+\Pi)^{-1}({\bf 1}-\Pi)({\bf 1}-\exp(v_{j}A\Delta
t))(\tilde{\mathbf{E}}\mathbf{{\mathcal{D}f}})_{j}^{n}.$

The last term can be split into two parts: the first one $-\Delta t\Delta
v(A+\Pi)^{-1}({\bf
1}-\Pi)\sum_{j}(\tilde{\mathbf{E}}\mathbf{{\mathcal{D}f}})_{j}^{n}$ vanishes
thanks to the conservative properties of the discrete operator
$\mathbf{{\mathcal{D}}}$ whereas the second one, we use (28) to get

$\displaystyle\Delta t\Delta v(A+\Pi)^{-1}\sum_{j}({\bf
1}-\Pi)\exp(v_{j}A\Delta
t)(\tilde{\mathbf{E}}\mathbf{{\mathcal{D}f}})_{j}^{n}\\!\\!$
$\displaystyle\\!\\!=\\!\\!$ $\displaystyle\\!\\!\Delta
v(A+\Pi)^{-1}\sum_{j}({\bf 1}-\Pi)(\exp(v_{j}A\Delta
t)\mathbf{f}_{j}^{n}-\mathbf{f}_{j}^{n+1})$ $\displaystyle=$
$\displaystyle\Delta v(A+\Pi)^{-1}\sum_{j}\Big{[}({\bf
1}-\Pi)\exp(v_{j}A\Delta t)\mathbf{f}_{j}^{n}-({\bf
1}-\Pi)\mathbf{f}_{j}^{n+1}\Big{]}.$

Finally, we then get

$\displaystyle\mathbf{E}^{n+1}$ $\displaystyle=$
$\displaystyle\mathbf{E}^{n}+\Delta v\sum_{j}(A+\Pi)^{-1}({\bf 1}-\Pi)({\bf
1}-\exp(v_{j}A\Delta t))\mathbf{f}_{j}^{n}$ (29) $\displaystyle+\Delta
v(A+\Pi)^{-1}\sum_{j}\Big{[}({\bf 1}-\Pi)\exp(v_{j}A\Delta
t)\mathbf{f}_{j}^{n}-({\bf 1}-\Pi)\mathbf{f}_{j}^{n+1}\Big{]}$
$\displaystyle=$ $\displaystyle\mathbf{E}^{n}+\Delta
v\sum_{j}(A+\Pi)^{-1}(({\bf 1}-\Pi)\mathbf{f}_{j}^{n}-({\bf
1}-\Pi)\mathbf{f}_{j}^{n+1}).$

By induction, if the discrete Poisson equation is satisfied at iteration $n$,
then it is satisfied at iteration $n+1$ and the proof is complete.

∎

######  Remark 4.

In this remark, we discuss how to compute the matrix $\Pi$. As mentioned in
Remark 2, when $(k+1)N_{x}$ is odd (referred as odd case), $0$ is a single
eigenvalue and the associated eigenvector $u_{1}$ corresponds to the constant
function in the DG space (see Appendix C). When $(k+1)N_{x}$ is even, $0$ has
a double multiplicity and one has to find another eigenvector $u_{2}$ (see
Appendix C). Once we get the eigenvectors $u_{1},\
u_{2}\in\mathbb{R}^{(k+1)N_{x}}$ of $A$ associated to the eigenvalue $0$, then
$\Pi$ is given by $\Pi x=\langle x,u_{1}\rangle u_{1}+\langle x,u_{2}\rangle
u_{2}=(u_{1}^{T}\otimes u_{1}+u_{2}^{T}\otimes u_{2})x$, for all
$x\in\mathbb{R}^{(k+1)N_{x}}$, with the Kronecker product $\otimes$ and
$(v^{T}\otimes v)_{i,j}=v_{i}v_{j}$. Some examples are given in Appendix C.

###  3.2 Vlasov Maxwell equations 1dx-2dv

In this part, We consider the following Vlasov-Maxwell 1dx-2dv model satisfied
by $f(t,x,v_{1},v_{2})$, $E_{1}(t,x)$, $E_{2}(t,x)$, $B(t,x)$, with $t\geq
0,x\in[0,L](L>0)$ and $(v_{1},v_{2})\in\mathbb{R}^{2}$

$\left\\{\begin{aligned}
&\partial_{t}f+v_{1}\partial_{x}f+E_{1}\partial_{v_{1}}f+E_{2}\partial_{v_{2}}f+B(v_{2}\partial_{v_{1}}f-v_{1}\partial_{v_{2}}f)=0,\\\
&\partial_{t}B=-\partial_{x}E_{2},\\\
&\partial_{t}E_{1}=-\int_{\mathbb{R}^{2}}v_{1}f\mathrm{d}v_{1}\mathrm{d}v_{2}+\bar{J}_{1},\\\
&\partial_{t}E_{2}=-\partial_{x}B-\int_{\mathbb{R}^{2}}v_{2}f\mathrm{d}v_{1}\mathrm{d}v_{2}+\bar{J}_{2},\end{aligned}\right.$
(30)

with the initial conditions $(f_{0}(x,v),E_{1}^{0}(x),E_{2}^{0}(x),B^{0}(x))$
such that the Poisson equation is satisfied initially
$\partial_{x}E_{1}^{0}(x)=\int_{\mathbb{R}^{2}}f_{0}(x,v)dv_{1}dv_{2}-\bar{\rho}$,
with
$\bar{\rho}=\frac{1}{L}\int_{0}^{L}\int_{\mathbb{R}^{2}}f_{0}(x,v)dxdv_{1}dv_{2}$
and periodic boundary conditions are imposed in space. Here
$\bar{J}_{i}=\frac{1}{L}\int_{0}^{L}\int_{\mathbb{R}^{2}}v_{i}fdxdv_{1}dv_{2}$
ensures that the electric fields are zero average in space.

####  3.2.1 Semi-discretization

We follow the lines of the above subsection and use a DG method in the space
direction $x$ and we consider a grid in the velocity direction
$v_{\ell,j_{\ell}}=-v_{\ell,\max}+j_{\ell}\Delta
v_{\ell},v_{\ell}\in[-v_{\ell,\max},v_{\ell,\max}],\ell=1,2,\Delta
v_{\ell}=2v_{\ell,\max}/{N_{v_{\ell}}},N_{v_{\ell}}\in\mathbb{N}^{\star}$. The
definitions of the different objects are a direct extension of the definitions
introduced in the previous part. Indeed, we denote by
$\mathbf{f}_{j_{1},j_{2}}\in\mathbb{R}^{(k+1)N_{x}}$ the DG coefficient vector
of $f_{j_{1},j_{2}}(t,x)=f_{h}(t,x,v_{j_{1}},v_{j_{2}})\approx
f(t,x,v_{j_{1}},v_{j_{2}})$ in space and evaluated at the velocity grid, and
$\mathbf{E_{1}},\mathbf{E_{2}},\mathbf{B}\in\mathbb{R}^{(k+1)N_{x}}$ are the
DG coefficient vectors of
$(E_{1,h},E_{2,h},B_{h})(t,x)\approx(E_{1},E_{2},B)(t,x)$. Moreover,
$(\tilde{\mathbfcal{F}}{\cal D}\mathbf{f})_{j_{1},j_{2}}$ (with
$\tilde{\mathbfcal{F}}_{j_{1},j_{2}}\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}(\mathbb{R})$)
is obtained as previously by a DG approximation of the nonlinear term $({\cal
F}\cdot\nabla_{v}f),\ \mbox{ with }{\cal
F}=({E_{1}}+{B}v_{2},{E_{2}}-{B}v_{1})$. For fixed indices $j_{1},j_{2}$, the
derivation of the numerical scheme is very similar to the 1dx-1dv case and we
then obtain the following semi-discretized (in space and velocity) scheme for
$j_{\ell}=1,\dots,N_{v_{\ell}}(\ell=1,2)$

$\left\\{\begin{aligned}
&\partial_{t}\mathbf{f}_{j_{1},j_{2}}=v_{1,j_{1}}A\mathbf{f}_{j_{1},j_{2}}-(\tilde{\mathbfcal{F}}{\cal
D}\mathbf{f})_{j_{1},j_{2}},\\\ &\partial_{t}\mathbf{B}=A\mathbf{E_{2}},\\\
&\partial_{t}\mathbf{E_{1}}=-\sum_{j_{1},j_{2}}v_{1,j_{1}}({\bf
1}-\Pi)\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2},\\\
&\partial_{t}\mathbf{E_{2}}=A\mathbf{B}-\sum_{j_{1},j_{2}}v_{2,j_{2}}({\bf
1}-\Pi)\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2},\end{aligned}\right.$
(31)

where $A$ is given by (8)-(9) and ${\bf 1}$ the identity matrix of size
$(k+1)N_{x}$ and $\Pi$ the projection matrix onto Ker$(A)$ introduced in the
previous part. Let denote
$U(t)=(\vec{\mathbf{f}}(t),\mathbf{B}(t),\mathbf{E_{1}}(t),\mathbf{E_{2}}(t))$,
where $\vec{\mathbf{f}}(t)\in\mathbb{R}^{(k+1)N_{x}N_{v_{1}}N_{v_{2}}}$
contains the DG coefficients
$\mathbf{f}_{j_{1},j_{2}}$,
$\mathbf{B}(t),\mathbf{E}_{1}(t),\mathbf{E}_{2}(t)\in\mathbb{R}^{(k+1)N_{x}}$
denote the DG coefficients of the electromagnetic fields. Using the above
notations and the following ones

$\displaystyle\vec{v_{\ell}}\in\mathbb{R}^{N_{v_{\ell}}},$
$\displaystyle(\vec{v_{\ell}})_{j}=v_{\ell,j}=-v_{\ell,\max}+j\Delta
v_{\ell},\;\;j=1,\dots,N_{v_{\ell}}\mbox{ and }\ell=1,2,$
$\displaystyle\vec{f}_{\star,j_{2}}\in\mathbb{R}^{(k+1)N_{x}N_{v_{1}}}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mbox{ for }j_{2}=1,\dots,N_{v_{2}}.$

We also introduce the following compact notations for the size of the
matrices: $\tilde{N}=(k+1)N_{x}$ and
$\tilde{N}_{1}=(k+1)N_{x}N_{v_{1}}=\tilde{N}N_{v_{1}}$. The system (31) can be
recast as

$\partial_{t}U=LU+N(U),$ (32)

with

$\displaystyle U$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{llll}\vec{\mathbf{f}}_{\star,1}\\\
\vec{\mathbf{f}}_{\star,2}\\\ \vdots\\\
\\!\\!\\!\vec{\mathbf{f}}_{\star,N_{v_{2}}}\\!\\!\\!\\!\\\ \mathbf{B}\\\
\mathbf{E_{2}}\\\
\mathbf{E_{1}}\end{array}\right),\;\;\;\;\;\;\;\;\;\;\;\;N(U)=\left(\begin{array}[]{llllccc}-(\tilde{\mathbfcal{F}}{\cal
D}\vec{\mathbf{f}})_{\star,1}\\\ -(\tilde{\mathbfcal{F}}{\cal
D}\vec{\mathbf{f}})_{\star,2}\\\ \vdots\\\ -(\tilde{\mathbfcal{F}}{\cal
D}\vec{\mathbf{f}})_{\star,N_{v_{2}}}\\\ {\bf 0}\\\ {\bf 0}\\\ {\bf
0}\end{array}\right),\;\;\;\;\;\;\;\;\;\;\;\;\vec{\mathbf{f}}_{\star,j_{2}}\in\mathbb{R}^{\tilde{N_{1}}},\forall
j_{2}=1,\dots,N_{v_{2}}$ (47) $\displaystyle L$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{lllllllccc}\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}\ldots&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&{\bf 0}_{\tilde{N_{1}},\tilde{N}}&{\bf
0}_{\tilde{N_{1}},\tilde{N}}&{\bf 0}_{\tilde{N_{1}},\tilde{N}}\\\ {\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}\ldots&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&\vdots&\vdots&\vdots\\\
\vdots&\ddots&\ddots&\vdots&\vdots&\vdots&\vdots\\\ {\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&\ldots&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}&{\bf
0}_{\tilde{N_{1}},\tilde{N}}&{\bf 0}_{\tilde{N_{1}},\tilde{N}}&{\bf
0}_{\tilde{N_{1}},\tilde{N}}\\\ {\bf 0}_{\tilde{N},\tilde{N_{1}}}&\ldots&{\bf
0}_{\tilde{N},\tilde{N_{1}}}&{\bf 0}_{\tilde{N},\tilde{N_{1}}}&{\bf
0}_{\tilde{N},\tilde{N}}&A&{\bf 0}_{\tilde{N},\tilde{N}}\\\ {\cal
E}_{2,1}&{\cal E}_{2,2}&\ldots&{\cal E}_{2,N_{v_{2}}}&A&{\bf
0}_{\tilde{N},\tilde{N}}&{\bf 0}_{\tilde{N},\tilde{N}}\\\ {\cal E}_{1}&{\cal
E}_{1}&\ldots&{\cal E}_{1}&{\bf 0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}&{\bf 0}_{\tilde{N},\tilde{N}}\\\ \end{array}\right)$
(55)

where diag$(\vec{v_{1}})\in\mathbb{M}_{N_{v_{1}},N_{v_{1}}}(\mathbb{R})$
denotes the diagonal matrix with $\vec{v_{1}}\in\mathbb{R}^{N_{v_{1}}}$ on its
diagonal, $\otimes$ denotes Kronecker product,
$\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}\in\mathbb{M}_{\tilde{N}_{1},\tilde{N}_{1}}$
and ${\bf 0}_{m,n}$ is the zero matrix with $m$ lines and $n$ columns.
Moreover, the matrices ${\cal E}_{1},{\cal
E}_{2,j_{2}}\in\mathbb{M}_{\tilde{N},\tilde{N}_{1}}(\mathbb{R})$ for
$j_{2}=1,\dots,N_{v_{2}}$ are defined by

$\displaystyle{\cal E}_{1}$ $\displaystyle=$ $\displaystyle-\Delta v_{1}\Delta
v_{2}\vec{v_{1}}\otimes({\bf 1}-\Pi),$ $\displaystyle{\cal E}_{2,j_{2}}$
$\displaystyle=$ $\displaystyle-\Delta v_{1}\Delta
v_{2}v_{2,j_{2}}\mathbbm{1}\otimes({\bf 1}-\Pi)\;\;\mbox{ with
}\;\;\mathbbm{1}=(1,\dots,1)\in\mathbb{R}^{N_{v_{1}}}.$

The size of the matrix $L$ is
$(\tilde{N_{1}}N_{v_{2}}+3\tilde{N})\times(\tilde{N_{1}}N_{v_{2}}+3\tilde{N})=(3+N_{v_{1}}N_{v_{2}})(k+1)N_{x}\times(3+N_{v_{1}}N_{v_{2}})(k+1)N_{x}$.
Even if it is a large matrix, one can see that $L$ is sparse which will help
to compute its exponential.

####  3.2.2 Time discretization

We now study the time discretization of (32) and as previously, we will design
an exponential scheme. To do so, the discrete unknown $U^{n}\approx
U(t^{n})\;(t^{n}=n\Delta t,\Delta t>0)$ is updated by

$U^{n+1}=\exp(\Delta tL)U^{n}+\Delta t\exp(\Delta tL)N(U^{n}),$ (56)

and one has to compute $\exp(\Delta tL)$. We give the following proposition to
show the representation of $\exp(\Delta tL)$.

######  Proposition 4.

The exponential of the matrix $L$ given by (55) is given by

$\exp(L\Delta t)=\left(\begin{array}[]{lllllllccc}e^{\Delta
t\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}}&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}\ldots&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&{\bf 0}_{\tilde{N_{1}},\tilde{N}}&{\bf
0}_{\tilde{N_{1}},\tilde{N}}&{\bf 0}_{\tilde{N_{1}},\tilde{N}}\\\ {\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&e^{\Delta
t\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}}&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}\ldots&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&\vdots&\vdots&\vdots\\\
\vdots&\ddots&\ddots&\vdots&\vdots&\vdots&\vdots\\\ {\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&\ldots&{\bf
0}_{\tilde{N_{1}},\tilde{N_{1}}}&e^{\Delta
t\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}}&{\bf
0}_{\tilde{N_{1}},\tilde{N}}&{\bf 0}_{\tilde{N_{1}},\tilde{N}}&{\bf
0}_{\tilde{N_{1}},\tilde{N}}\\\ {}_{e}{\cal B}_{1}&{}_{e}{\cal
B}_{2}&\ldots&{}_{e}{\cal B}_{N_{v_{2}}}&\frac{e^{A\Delta t}+e^{-A\Delta
t}}{2}&\frac{e^{A\Delta t}-e^{-A\Delta t}}{2}&{\bf 0}_{\tilde{N},\tilde{N}}\\\
{}_{e}{\cal E}_{2,1}&{}_{e}{\cal E}_{2,2}&\ldots&{}_{e}{\cal
E}_{2,N_{v_{2}}}&\frac{e^{A\Delta t}-e^{-A\Delta t}}{2}&\frac{e^{A\Delta
t}+e^{-A\Delta t}}{2}&{\bf 0}_{\tilde{N},\tilde{N}}\\\ {}_{e}{\cal
E}_{1}&{}_{e}{\cal E}_{1}&\ldots&{}_{e}{\cal E}_{1}&{\bf
0}_{\tilde{N},\tilde{N}}&{\bf 0}_{\tilde{N},\tilde{N}}&{\bf
1}\end{array}\right)$

where ${}_{e}{\cal B}_{j_{2}},{}_{e}{\cal E}_{2,j_{2}},{}_{e}{\cal
E}_{1}\in\mathbb{M}_{\tilde{N},\tilde{N}_{1}}(\mathbb{R})$ for
$j_{2}=1,\dots,N_{v_{2}}$ are given by

$\displaystyle{}_{e}{\cal B}_{j_{2}}$ $\displaystyle=$ $\displaystyle\Delta
v_{1}\Delta v_{2}\;v_{2,j_{2}}(A+\Pi)^{-1}\vec{\alpha},$
$\displaystyle{}_{e}{\cal E}_{2,j_{2}}$ $\displaystyle=$ $\displaystyle\Delta
v_{1}\Delta v_{2}\;v_{2,j_{2}}(A+\Pi)^{-1}\vec{\beta},$
$\displaystyle{}_{e}{\cal E}_{1}$ $\displaystyle=$ $\displaystyle\Delta
v_{1}\Delta v_{2}(A+\Pi)^{-1}(\mathbbm{1}\otimes{\bf 1}-e^{\Delta
t\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}})$

where the matrices $\vec{\alpha},\vec{\beta}$ are given by

$\displaystyle\vec{\alpha}$ $\displaystyle=$
$\displaystyle[\alpha_{1},\alpha_{2},\dots,\alpha_{N_{v_{1}}}]\in\mathbb{M}_{\tilde{N},\tilde{N_{1}}}(\mathbb{R})\mbox{
with }\alpha_{j_{1}}=\Big{[}\frac{-e^{A\Delta
t}}{2(1-v_{1,j_{1}})}-\frac{e^{-A\Delta t}}{2(1+v_{1,j_{1}})}+\frac{e^{A\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})}\Big{]},$ $\displaystyle\vec{\beta}$
$\displaystyle=$
$\displaystyle[\beta_{1},\beta_{2},\dots,\beta_{N_{v_{1}}}]\in\mathbb{M}_{\tilde{N},\tilde{N_{1}}}(\mathbb{R})\mbox{
with }\beta_{j_{1}}=\Big{[}\frac{-e^{A\Delta
t}}{2(1-v_{1,j_{1}})}+\frac{e^{-A\Delta
t}}{2(1+v_{1,j_{1}})}+\frac{v_{1,j_{1}}e^{A\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})}\Big{]}.$

###### Proof.

First the $(\tilde{N_{1}}N_{v_{2}})\times(\tilde{N_{1}}N_{v_{2}})$ block is
diagonal and the diagonal part is $e^{\Delta
t\small{\mbox{diag}(\vec{v_{1}})\\!\otimes\\!A}}$ ($N_{v_{2}}$ times). Second,
the $3\times 3$ right bottom block corresponds to the homogeneous Maxwell
equations. Its exponential can be computed and is equal to

$\exp\left(\begin{array}[]{lll}{\bf 0}_{\tilde{N},\tilde{N}}&A\Delta t&{\bf
0}_{\tilde{N},\tilde{N}}\\\ A\Delta t&{\bf 0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}\\\ {\bf 0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}\end{array}\right)=\left(\begin{array}[]{lll}\frac{\exp(A\Delta
t)+\exp(-A\Delta t)}{2}&\frac{\exp(A\Delta t)-\exp(-A\Delta t)}{2}&{\bf
0}_{\tilde{N},\tilde{N}}\\\ \frac{\exp(A\Delta t)-\exp(-A\Delta
t)}{2}&\frac{\exp(A\Delta t)+\exp(-A\Delta t)}{2}&{\bf
0}_{\tilde{N},\tilde{N}}\\\ {\bf 0}_{\tilde{N},\tilde{N}}&{\bf
0}_{\tilde{N},\tilde{N}}&{\bf 1}\end{array}\right).$

Finally, we compute the three last block lines of $\exp(L\Delta t)$.

##### Computation of ${}_{e}{\cal E}_{1}$ : solve $\mathbf{E_{1}}$

First, we have for $\mathbf{f}_{j_{1},j_{2}}(t)$

$\mathbf{f}_{j_{1},j_{2}}(t)=e^{v_{1,j_{1}}A(t-t^{n})}\mathbf{f}_{j_{1},j_{2}}(t^{n}),$

which enables to compute $\mathbf{E_{1}}(t^{n+1})$

$\displaystyle\mathbf{E_{1}}(t^{n+1})$ $\displaystyle=$
$\displaystyle\mathbf{E_{1}}(t^{n})-\sum_{j_{1},j_{2}}\int_{t^{n}}^{t^{n+1}}({\bf
1}-\Pi)e^{v_{1,j_{1}}A(t-t^{n})}dtv_{1,j_{1}}\mathbf{f}_{j_{1},j_{2}}(t^{n})\Delta
v_{1}\Delta v_{2}$ $\displaystyle=$ $\displaystyle\mathbf{E_{1}}(t^{n})+\Delta
v_{1}\Delta v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-e^{v_{1,j_{1}}A\Delta
t})\mathbf{f}_{j_{1},j_{2}}(t^{n}),$

from which we can thus deduce the last line of the exponential of the matrix.

##### Computation of ${}_{e}{\cal B}_{j_{2}},{}_{e}{\cal E}_{2,j_{2}}$: solve
$\mathbf{B},\mathbf{E_{2}}$

Next, we focus on the calculation of $\mathbf{E_{2}}(t^{n+1})$ and
$\mathbf{B}(t^{n+1})$ from known initial conditions $\mathbf{E_{2}}(t^{n})$
and $\mathbf{B}(t^{n})$. Let write down the equations for
$(\mathbf{E_{2}},\mathbf{B})(t)$ with $t\in[t^{n},t^{n+1}]$

$\displaystyle\frac{d}{dt}\mathbf{E_{2}}(t)$ $\displaystyle=$ $\displaystyle
A\mathbf{B}(t)-\sum_{j_{1},j_{2}}({\bf
1}-\Pi)e^{v_{1,j_{1}}A(t-t^{n})}v_{2,j_{1}}\mathbf{f}_{j_{1},j_{2}}(t^{n})\Delta
v_{1}\Delta v_{2}$ $\displaystyle\frac{d}{dt}\mathbf{B}(t)$ $\displaystyle=$
$\displaystyle A\mathbf{E_{2}}(t)$

which can be rewritten as $\frac{dU}{dt}=MU+R$ with
$U(t)=(\mathbf{E_{2}}(t),\mathbf{B}(t))$ and

$\displaystyle M$ $\displaystyle=$ $\displaystyle\left(\begin{array}[]{ll}{\bf
0}_{\tilde{N},\tilde{N}}&A\\\ A&{\bf
0}_{\tilde{N},\tilde{N}}\end{array}\right),\;\;R(t)=\left(\begin{array}[]{ll}R_{1}(t)\\\
{\bf 0}_{\tilde{N},1}\end{array}\right),$ (61) $\displaystyle R_{1}(t)$
$\displaystyle=$ $\displaystyle-\sum_{j_{1},j_{2}}({\bf
1}-\Pi)e^{v_{1,j_{1}}A(t-t^{n})}v_{2,j_{1}}\mathbf{f}_{j_{1},j_{2}}(t^{n})\Delta
v_{1}\Delta v_{2}$ (62)

Thus, one can write the variation of constant formula

$U(t^{n+1})=e^{M\Delta
t}U(t^{n})+\int_{t^{n}}^{t^{n+1}}e^{-M(t-t^{n+1})}R(t)dt.$ (63)

First, $e^{M\Delta t}$ reads as, using its definition (62)

$e^{M\Delta t}=\frac{1}{2}\left(\begin{array}[]{ll}{e^{A\Delta t}+e^{-A\Delta
t}}&{e^{A\Delta t}-e^{-A\Delta t}}\\\ {e^{A\Delta t}-e^{-A\Delta
t}}&{e^{A\Delta t}+e^{-A\Delta t}}\end{array}\right).$

Second, one has to compute the integral term in (63)

$\displaystyle\int_{t^{n}}^{t^{n+1}}e^{-M(t-t^{n+1})}R(t)dt$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\begin{array}[]{ll}\int_{t^{n}}^{t^{n+1}}\Big{[}e^{-A(t-t^{n+1})}+e^{A(t-t^{n+1})}\Big{]}R_{1}(t)dt\\\
\int_{t^{n}}^{t^{n+1}}\Big{[}e^{-A(t-t^{n+1})}-e^{A(t-t^{n+1})}\Big{]}R_{1}(t)dt\\\
\end{array}\right)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(\begin{array}[]{ll}\int_{t^{n}}^{t^{n+1}}\sum_{j_{1},j_{2}}\Big{[}e^{-A(t-t^{n+1})}+e^{A(t-t^{n+1})}\Big{]}e^{v_{1,j_{1}}A(t-t^{n})}v_{2,j_{1}}{f}^{n}_{k,j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}dt\\\
\int_{t^{n}}^{t^{n+1}}\sum_{j_{1},j_{2}}\Big{[}e^{-A(t-t^{n+1})}-e^{A(t-t^{n+1})}\Big{]}e^{v_{1,j_{1}}A(t-t^{n})}v_{2,j_{1}}{f}^{n}_{k,j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}dt\end{array}\right)$ $\displaystyle=$
$\displaystyle-\left(\begin{array}[]{ll}\sum_{j_{1},j_{2}}{\cal
I}_{1}\;v_{2,j_{1}}{f}^{n}_{k,j_{1},j_{2}}\Delta v_{1}\Delta v_{2}\\\
\sum_{j_{1},j_{2}}{\cal I}_{2}\;v_{2,j_{1}}{f}^{n}_{k,j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}\end{array}\right)$

where ${\cal I}_{1},{\cal I}_{2}$ are given by

$\displaystyle{\cal I}_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t^{n}}^{t^{n+1}}\Big{[}e^{-A(t-t^{n+1})}+e^{A(t-t^{n+1})}\Big{]}e^{v_{1,j_{1}}A(t-t^{n})}dt$
$\displaystyle=$ $\displaystyle\frac{(A+\Pi)^{-1}e^{A\Delta
t}}{2(1-v_{1,j_{1}})}-\frac{(A+\Pi)^{-1}e^{-A\Delta
t}}{2(1+v_{1,j_{1}})}-\frac{(A+\Pi)^{-1}v_{1,j_{1}}e^{A\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})},$ $\displaystyle{\cal I}_{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{t^{n}}^{t^{n+1}}\Big{[}e^{-A(t-t^{n+1})}-e^{A(t-t^{n+1})}\Big{]}e^{v_{1,j_{1}}A(t-t^{n})}dt$
$\displaystyle=$ $\displaystyle\frac{(A+\Pi)^{-1}e^{A\Delta
t}}{2(1-v_{1,j_{1}})}+\frac{(A+\Pi)^{-1}e^{-A\Delta
t}}{2(1+v_{1,j_{1}})}-\frac{(A+\Pi)^{-1}e^{A\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})}.$

Inserting these calculations in (63) leads to the following expression for
$\mathbf{E_{2}}(t^{n+1})$ and $\mathbf{B}(t^{n+1})$

$\displaystyle\mathbf{E_{2}}(t^{n+1})$ $\displaystyle=$
$\displaystyle\frac{1}{2}(e^{A\Delta t}+e^{-A\Delta
t})\mathbf{E_{2}}(t^{n})+\frac{1}{2}(e^{A\Delta t}-e^{-A\Delta
t})\mathbf{B}(t^{n})+\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}v_{2,j_{2}}\beta_{j_{1}}\mathbf{f}_{j_{1},j_{2}}(t^{n})$
$\displaystyle\mathbf{B}(t^{n+1})$ $\displaystyle=$
$\displaystyle\frac{1}{2}(e^{A\Delta t}-e^{-A\Delta
t})\mathbf{E_{2}}(t^{n})+\frac{1}{2}(e^{A\Delta t}+e^{-A\Delta
t})\mathbf{B}(t^{n})+\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}v_{2,j_{2}}\alpha_{j_{1}}\mathbf{f}_{j_{1},j_{2}}(t^{n})$

where
$\vec{\beta}=[\beta_{1},\beta_{2},\dots,\beta_{N_{v_{1}}}]\in\mathbb{M}_{\tilde{N}\times\tilde{N_{1}}}(\mathbb{R})$
and
$\vec{\alpha}=[\alpha_{1},\alpha_{2},\dots,\alpha_{N_{v_{1}}}]\in\mathbb{M}_{\tilde{N}\times\tilde{N_{1}}}(\mathbb{R})$
are given by

$\displaystyle\beta_{j_{1}}$ $\displaystyle=$ $\displaystyle\frac{-e^{A\Delta
t}}{2(1-v_{1,j_{1}})}+\frac{e^{-A\Delta
t}}{2(1+v_{1,j_{1}})}+\frac{v_{1,j_{1}}e^{A\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})}$ $\displaystyle\alpha_{j_{1}}$
$\displaystyle=$ $\displaystyle\frac{-e^{A\Delta
t}}{2(1-v_{1,j_{1}})}-\frac{e^{-A\Delta t}}{2(1+v_{1,j_{1}})}+\frac{e^{A\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})},$

which conclude the proof.

∎

Shared the same spirit with Vlasov-Ampère equation, we have the following
discretized Poisson equation and error estimate for Vlasov-Maxwell equation
with exponential Lawson RK DG FD discretization.

######  Proposition 5.

The exponential DG method (56) where the exponential is given in Prop 4 (and
its generalization to high order Lawson Runge-Kutta) satisfied by
$U^{n}=(\vec{\bf{f}},{\bf B},{\bf E}_{2},{\bf E}_{1})^{n}$ preserves the
following Poisson equation

$(A+\Pi)\mathbf{E}_{1}^{n}=-\sum_{j_{1},j_{2}}({\bf
1}-\Pi)\mathbf{f}^{n}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2},\;\;\;\forall
n\in\mathbb{N}^{\star},$

provided that it is satisfied at the initial time $n=0$. Here, $A$ is the DG
matrix given by (8)-(9), $\Pi$ the orthogonal projection onto Ker$(A)$ and
${\bf{1}}$ is the identity matrix of size $\tilde{N}=(k+1)N_{x}$.

###### Proof.

We present the proof for first order Lawson case (forward Euler), the proof
can be generalized to arbitrary explicit Runge-Kutta scheme. First, we assume
the Poisson equation $(A+\Pi)\mathbf{E}_{1}^{0}=-\sum_{j_{1},j_{2}}({\bf
1}-\Pi)\mathbf{f}^{0}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2}$ holds at the
initial time.

Next, from the scheme (56) and Prop 4, we have

$\mathbf{f}_{j_{1},j_{2}}^{n+1}=e^{v_{1,j_{1}}A\Delta
t}\mathbf{f}_{j_{1},j_{2}}^{n}-\Delta te^{v_{1,j_{1}}A\Delta
t}(\tilde{\mathbfcal{F}}{\cal D}\mathbf{f})_{j_{1},j_{2}}^{n}.$

Regarding the $\mathbf{E}_{1}$ component, we have

$\displaystyle\mathbf{E}^{n+1}_{1}$ $\displaystyle=\mathbf{E}^{n}_{1}+\Delta
v_{1}\Delta v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-e^{v_{1,j_{1}}A\Delta
t})\mathbf{f}_{j_{1},j_{2}}^{n}$ $\displaystyle-\Delta t\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-e^{v_{1,j_{1}}A\Delta
t})(\tilde{\mathbfcal{F}}{\cal D}\mathbf{f})_{j_{1},j_{2}}^{n}$
$\displaystyle=\mathbf{E}^{n}_{1}+\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-\Pi)({\bf 1}-e^{v_{1,j_{1}}A\Delta
t})\mathbf{f}_{j_{1},j_{2}}^{n}$ $\displaystyle-\Delta t\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-\Pi)({\bf 1}-e^{v_{1,j_{1}}A\Delta
t})(\tilde{\mathbfcal{F}}{\cal D}\mathbf{f})_{j_{1},j_{2}}^{n}$
$\displaystyle=\mathbf{E}^{n}_{1}+\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-\Pi)({\bf 1}-e^{v_{1,j_{1}}A\Delta
t})\mathbf{f}_{j_{1},j_{2}}^{n}$ $\displaystyle-\Delta t\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}({\bf 1}-\Pi)\sum_{j_{1},j_{2}}(\tilde{\mathbfcal{F}}{\cal
D}\mathbf{f})_{j_{1},j_{2}}^{n}$ $\displaystyle+\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}({\bf 1}-\Pi)(e^{v_{1,j_{1}}A\Delta
t}\mathbf{f}_{j_{1},j_{2}}^{n}-\mathbf{f}_{j_{1},j_{2}}^{n+1})$
$\displaystyle=\mathbf{E}^{n}_{1}+\Delta v_{1}\Delta
v_{2}(A+\Pi)^{-1}\sum_{j_{1},j_{2}}(({\bf
1}-\Pi)\mathbf{f}_{j_{1},j_{2}}^{n}-({\bf
1}-\Pi)\mathbf{f}_{j_{1},j_{2}}^{n+1}).$

By induction, if the discrete Poisson equation is satisfied at iteration $n$,
then it is satisfied at iteration $n+1$ and the proof is complete.

∎

##  4 Vlasov-Maxwell 2dx-2dv

We finally consider the 2dx-2dv Vlasov-Maxwell model satisfied by
$f(t,x,y,v_{1},v_{2}),E_{1}(t,x,y),E_{2}(t,x,y),B(t,x,y)$, with $t\geq
0,(x,y)\in[0,L_{x}]\times[0,L_{y}]$ and $(v_{1},v_{2})\in\mathbb{R}^{2}$

$\left\\{\begin{aligned}
&\partial_{t}f+v_{1}\partial_{x}f+v_{2}\partial_{y}f+E_{1}\partial_{v_{1}}f+E_{2}\partial_{v_{2}}f+B(v_{2}\partial_{v_{1}}f-v_{1}\partial_{v_{2}}f)=0,\\\
&\partial_{t}B=\partial_{y}E_{1}-\partial_{x}E_{2},\\\
&\partial_{t}E_{1}=\partial_{y}B-\int_{\mathbb{R}^{2}}v_{1}f\mathrm{d}v_{1}\mathrm{d}v_{2}+\bar{J}_{1},\\\
&\partial_{t}E_{2}=-\partial_{x}B-\int_{\mathbb{R}^{2}}v_{2}f\mathrm{d}v_{1}\mathrm{d}v_{2}+\bar{J}_{2},\\\
&\partial_{x}E_{1}+\partial_{y}E_{2}=\int_{\mathbb{R}^{2}}f\mathrm{d}v_{1}\mathrm{d}v_{2}-\bar{\rho},\;\;\;\partial_{x}B+\partial_{y}B=0,\end{aligned}\right.$
(67)

with initial conditions
$(f_{0}(x,y,v_{1},v_{2}),E^{0}_{1}(x,y),E_{2}^{0}(x,y),B(x,y))$ such that the
Poisson equation is satisfied
$\partial_{x}E_{1}^{0}+\partial_{y}E_{2}^{0}=\int_{\mathbb{R}^{2}}f_{0}\mathrm{d}v_{1}\mathrm{d}v_{2}-\bar{\rho}$
and periodic boundary conditions are considered in space. here,
$\bar{J}_{i}=\frac{1}{L_{x}L_{y}}\int_{L_{x}\times
L_{y}}\int_{\mathbb{R}^{2}}v_{i}fdxdydv_{1}dv_{2}$,
$\bar{\rho}=\frac{1}{L_{x}L_{y}}\int_{L_{x}\times
L_{y}}\int_{\mathbb{R}^{2}}fdxdydv_{1}dv_{2}$,
$E=(E_{1},E_{2}),B=(B_{1},B_{2}),\nabla=(\partial_{x},\partial_{y})$.

###  4.1 Exponential DG discretization

Here we apply 2D DG method in $(x,y)$ direction and consider the
discretization on Cartesian meshes with a rectangular triangulation
$I_{i}\times I_{j}$. We define the space $V_{h}$ as the space of tensor
product piece-wise polynomials of degree at most $k$ in each variable on every
element, i.e. $V_{h}^{k}=\\{v_{h}:v_{h}|_{I_{i}\times I_{j}}\in
Q^{k}(I_{i}\times I_{j})\\}$, where $Q^{k}(I_{i}\times I_{j})$ is the space of
tensor products of one dimensional polynomials of degree up to $k$.

We follow the lines of the above subsections: we use a DG method in the 2D
space direction $(x,y)$ (with $N_{x}$ (resp. $N_{y}$) cells in direction $x$
(resp. $y$) and a grid in the velocity direction
$v_{\ell,j_{\ell}}=v_{\ell,\min}+j_{\ell}\Delta
v_{\ell},\ell=1,2,\;\;j_{\ell}=1,\dots,N_{v}$. The 2D DG approximation for $f$
is represented as (for $j_{1},j_{2}=1,\dots,N_{v}$)

$f(t,x,y,v_{j_{1}},v_{j_{2}})\approx
f_{h}(t,x,y,v_{j_{1}},v_{j_{2}})=\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\sum_{m=0}^{k}\sum_{n=0}^{k}{f}_{ij}^{mn}(t,v_{j_{1}},v_{j_{2}})\xi_{i}^{m}(x)\xi_{j}^{n}(y).$

For simplicity, we only present the 2D DG discretization for linear part of
equation of $f$ obtained by multiplying the Vlasov equation by
$\xi^{\ell}(x)\xi^{s}(y)$ (for $\ell=0,\dots,k$ and $s=0,\dots,k$) on
$(x,y)\in I_{i}\times I_{j}$ (for $i=1,2,...N_{x}$ and $j=1,2,...N_{y}$):

$\displaystyle\sum_{m=0}^{k}\sum_{n=0}^{k}\left[\partial_{t}{f}_{ij}^{mn}(t,v_{j_{1}},v_{j_{2}})(\xi^{m},\xi^{\ell})_{I_{i}}(\xi^{n},\xi^{s})_{I_{j}}\right]$
(68)
$\displaystyle-\sum_{m=0}^{k}\sum_{n=0}^{k}\left[v_{j_{1}}{f}_{ij}^{mn}(t,v_{j_{1}},v_{j_{2}})(\xi^{m},\partial_{x}\xi^{\ell})_{I_{i}}(\xi^{n},\xi^{s})_{I_{j}}+v_{j_{2}}{f}_{ij}^{mn}(t,v_{j_{1}},v_{j_{2}})(\xi^{m},\xi^{\ell})_{I_{i}}(\xi^{n},\partial_{y}\xi^{s})_{I_{j}}\right]$
$\displaystyle+\sum_{m=0}^{k}\sum_{n=0}^{k}\left[(v_{j_{1}}\Big{[}\\{f_{h}(t,x,y,v_{j_{1}},v_{j_{2}}\\}\xi^{\ell}\Big{]}^{i+\frac{1}{2}}_{i-\frac{1}{2}},\xi^{s})_{I_{j}}+(v_{j_{2}}\Big{[}\\{f_{h}(t,x,y,v_{j_{1}},v_{j_{2}}\\}\xi^{s}\Big{]}^{j+\frac{1}{2}}_{j-\frac{1}{2}},\xi^{\ell})_{I_{i}}\right]=0,$

where we used the central fluxes in $x$ and $y$, that is for the $x$ direction
$\\{f_{h}(t,x,y,v_{j_{1}},v_{j_{2}})\\}|_{x_{i\pm
1/2}}=\frac{1}{2}(f_{h}(t,x_{i\pm
1/2}^{+},y,v_{j_{1}},v_{j_{2}})+f_{h}(t,x_{i\pm
1/2}^{-},y,v_{j_{1}},v_{j_{2}}))$.

We consider $\mathbf{f}_{j_{1},j_{2}}\in\mathbb{R}^{(k+1)^{2}N_{x}N_{y}}$ the
vector containing the degree of freedom $f^{m,n}_{i,j}$

$\mathbf{f}_{j_{1},j_{2}}=[f^{0,0}_{1,1},f^{1,0}_{1,1},\dots,f^{k,0}_{1,1},\dots,f^{k,0}_{N_{x},1},f^{0,1}_{1,1},\dots,f^{k,k}_{N_{x},N_{y}}]_{j_{1},j_{2}}^{T}$
(69)

and
$\mathbf{E_{1}},\mathbf{E_{2}},\mathbf{B}\in\mathbb{R}^{(k+1)^{2}N_{x}N_{y}}$
the vectors (defined as (69)) containing the DG degree of freedom of the
electromagnetic fields $(E_{1},E_{2},B)$. Finally, we introduce
$(\tilde{\mathbfcal F}{\cal D}\mathbf{f})_{j_{1},j_{2}}$ the DG approximation
of the nonlinear term $({\cal F}\cdot\nabla_{v}f)(v_{j_{1}},v_{j_{2}}),\
\mbox{ with }{\cal F}=({E_{1}}+{B}v_{2},{E_{2}}-{B}v_{1})$ using similar
techniques as in the 1dx case. With these notations, we have the following
semi-discretized scheme

$\left\\{\begin{aligned} &\partial_{t}\mathbf{f}_{j_{1},j_{2}}=({\bf
1}_{y}\otimes
v_{1,j_{1}}A^{x})\mathbf{f}_{j_{1},j_{2}}+(v_{2,j_{2}}A^{y}\otimes{\bf
1}_{x})\mathbf{f}_{j_{1},j_{2}}-(\tilde{\mathbfcal F}{\cal
D}\mathbf{f})_{j_{1},j_{2}}\\\ &\partial_{t}\mathbf{B}=-(A^{y}\otimes{\bf
1}_{x})\mathbf{E_{1}}+({\bf 1}_{y}\otimes A^{x})\mathbf{E_{2}},\\\
&\partial_{t}\mathbf{E_{1}}=-(A^{y}\otimes{\bf
1}_{x})\mathbf{B}-\sum_{j_{1},j_{2}}v_{1,j_{1}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2},\\\
&\partial_{t}\mathbf{E_{2}}=({\bf 1}_{y}\otimes
A^{x})\mathbf{B}-\sum_{j_{1},j_{2}}v_{2,j_{2}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2},\end{aligned}\right.$ (70)

where $A^{x}\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{x}}(\mathbb{R})$ and
$A^{y}\in\mathbb{M}_{(k+1)N_{y},(k+1)N_{y}}(\mathbb{R})$ are the matrices
coming from the DG semi-discretization in each space direction as before,
$\otimes$ is the Kronecker product, ${\bf P}=({\bf 1}_{y}-\Pi_{y})\otimes({\bf
1}_{x}-\Pi_{x})$ with ${\bf 1}_{x}$ (resp. ${\bf 1}_{y})$ the identity matrix
of size $(k+1)N_{x}$ (resp. $(k+1)N_{y}$) and $\Pi_{x}$ (resp. $\Pi_{y}$) the
projection onto Ker$(A^{x})$ (resp. Ker$(A^{y})$).

Before discussing the time discretization, we prove the following proposition.

######  Proposition 6.

The semi-discretized system (70) satisfied by $({{\bf f}}_{j_{1},j_{2}},{\bf
B},{\bf E}_{1},{\bf E}_{2})(t)$ preserves the following discretized Poisson
equation

$({\bf 1}_{y}\otimes A^{x})\mathbf{E_{1}}(t)+(A^{y}\otimes{\bf
1}_{x})\mathbf{E_{2}}(t)=-\sum_{j_{1},j_{2}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}(t)\Delta v_{1}\Delta v_{2},\;\;\;\forall t>0,$

provided that it is satisfied at time $t=0$.

###### Proof.

Let derive with respect to time the left hand side of the equality to get

$\displaystyle\partial_{t}(({\bf 1}_{y}\otimes
A^{x})\mathbf{E_{1}}+(A^{y}\otimes{\bf 1}_{x})\mathbf{E_{2}})$
$\displaystyle=$ $\displaystyle({\bf 1}_{y}\otimes
A^{x})\partial_{t}\mathbf{E_{1}}+(A^{y}\otimes{\bf
1}_{x})\partial_{t}\mathbf{E_{2}}$ $\displaystyle=$ $\displaystyle({\bf
1}_{y}\otimes A^{x})(-(A^{y}\otimes{\bf
1}_{x})\mathbf{B}-\sum_{j_{1},j_{2}}v_{1,j_{1}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2})$
$\displaystyle+(A^{y}\otimes{\bf 1}_{x})(({\bf 1}_{y}\otimes
A^{x})\mathbf{B}-\sum_{j_{1},j_{2}}v_{2,j_{2}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2})$ $\displaystyle=$
$\displaystyle\Big{[}(A^{y}\otimes{\bf 1}_{x})({\bf 1}_{y}\otimes A^{x})-({\bf
1}_{y}\otimes A^{x})(A^{y}\otimes{\bf 1}_{x})\Big{]}\mathbf{B}$
$\displaystyle-({\bf 1}_{y}\otimes A^{x})\sum_{j_{1},j_{2}}v_{1,j_{1}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2})$
$\displaystyle-(A^{y}\otimes{\bf 1}_{x})\sum_{j_{1},j_{2}}v_{2,j_{2}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2})$ $\displaystyle=$
$\displaystyle-\sum_{j_{1},j_{2}}\Big{(}\partial_{t}({\bf
P}\mathbf{f}_{j_{1},j_{2}})+{\bf P}(\tilde{\mathbfcal F}{\cal
D}\mathbf{f})_{j_{1},j_{2}}\Big{)}\Delta v_{1}\Delta v_{2}$ $\displaystyle=$
$\displaystyle-\partial_{t}\Big{(}\sum_{j_{1},j_{2}}{\bf
P}\mathbf{f}_{j_{1},j_{2}}\Delta v_{1}\Delta v_{2}\Big{)},$

where we used the identities

$\displaystyle(A^{y}\otimes{\bf 1}_{x})({\bf 1}_{y}\otimes A^{x})$
$\displaystyle=$ $\displaystyle(A^{y}{\bf 1}_{y})\otimes({\bf
1}_{x}A^{x})=A^{y}\otimes A^{x}$ $\displaystyle({\bf 1}_{y}\otimes
A^{x})(A^{y}\otimes{\bf 1}_{x})$ $\displaystyle=$ $\displaystyle({\bf
1}_{y}A^{y})\otimes(A^{x}{\bf 1}_{x})=A^{y}\otimes A^{x}$

to pass from the third to the fourth equality. Integrating in time the
obtained equality and assuming the discrete Poisson equation holds at time
$t=0$ leads to the result. ∎

We end this part by giving some elements on the time discretization. First, in
this case, it is difficult to compute the exponential of the linear part.
However, we can consider the exponential of the ${\bf f}_{j_{1},j_{2}}$ linear
part (which corresponds to the $(x,y)$ transport). Indeed, we observe from
(70) that this linear part writes

$\partial_{t}{\bf f}_{j_{1},j_{2}}=\Big{[}({\bf 1}_{y}\otimes
v_{1,j_{1}}A^{x})+(v_{2,j_{2}}A^{y}\otimes{\bf
1}_{x})\Big{]}\mathbf{f}_{j_{1},j_{2}}=\Big{[}v_{2,j_{2}}A^{y}\oplus
v_{1,j_{1}}A^{x}\Big{]}{\bf f}_{j_{1},j_{2}},$

where we used the definition of the Kronecker sum $\oplus$. The exact solution
can be then written as

${\bf f}_{j_{1},j_{2}}(t)=\exp\Big{(}(v_{2,j_{2}}A^{y}\oplus
v_{1,j_{1}}A^{x})\,t\Big{)}{\bf
f}_{j_{1},j_{2}}(0)=\exp\Big{(}v_{2,j_{2}}tA^{y}\oplus
v_{1,j_{1}}tA^{x}\Big{)}{\bf f}_{j_{1},j_{2}}(0).$

It is well known that the exponential of a matrix with Kronecker sum structure
is equal to the Kronecker product of the exponentials that is

${\bf
f}_{j_{1},j_{2}}(t)=\exp(v_{2,j_{2}}tA^{y})\otimes\exp(v_{1,j_{1}}tA^{x}){\bf
f}_{j_{1},j_{2}}(0),$

which can be recast using the vec operation as

$\Large{\textgoth{f}}_{j_{1},j_{2}}(t)=\exp(v_{2,j_{2}}tA^{y})\Large{\textgoth{f}}_{j_{1},j_{2}}(0)\exp(\Delta
v_{1,j_{1}}t(A^{x})^{T}),$ (71)

where ${\bf f}_{j_{1},j_{2}}(s)$=vec$(\Large{\textgoth{f}}_{j_{1},j_{2}}(s))$
(for $s=0,t$ denotes the vectorization operation which takes the matrix
$\Large{\textgoth{f}}_{j_{1},j_{2}}\in\mathbb{M}_{(k+1)N_{x},(k+1)N_{y}}(\mathbb{R})$
as entry and gives the vector ${\bf
f}_{j_{1},j_{2}}\in\mathbb{R}^{(k+1)^{2}N_{x}N_{y}}$ as a result. This means
that the update of $f$ requires matrix-vector products operations that only
involves to assembly exponential of matrices $A^{x}$ and $A^{y}$ which are
computed from the one-dimensional case (see (8)). Moreover, these matrix-
vector products calculations can be performed in a very efficient way. This
nice property has been exploited in the literature to design efficient
routines for computing matrix exponentials [6, 13, 24].

######  Remark 5.

The semi-discretized Vlasov-Maxwell system (70) can be degenerated to a semi-
discretization of the Vlasov-Poisson system satisfied by $(f,E_{1},E_{2})$. In
this case, the Lawson scheme only applies to the unknown ${\bf f}$ and then
requires the calculation of $\exp(v_{2,j_{2}}tA^{y}\oplus v_{1,j_{1}}tA^{x})$
which can be performed efficiently thanks to (71). The update of the electric
field is performed using the Poisson equation thanks to the updated ${\bf f}$.

###  4.2 Fourier based space discretization

In this part, we consider Fourier in space combined with finite differences in
velocity to semi-discretize the Vlasov-Maxwell system (67) and we will see
that in this case, it will be possible to compute explicitely the exponential
of the linear part.

Denoting $\hat{f}_{k_{x},k_{y},j_{1},j_{2}}$ the Fourier coefficient of $f$ in
space and evaluated at the velocity grid introduced previously ($k_{x},k_{y}$
being the Fourier variables),
$\hat{E}_{1,k_{x},k_{y}},\hat{E}_{2,k_{x},k_{y}},\hat{B}_{k_{x},k_{y}}$ the
Fourier coefficients of the electromagnetic fields $({E}_{1},{E}_{2},{B})$,
and introducing the force term ${\cal F}=(E_{1}+Bv_{2},E_{2}-Bv_{1})$, we get
the following semi-discretized scheme

$\left\\{\begin{aligned}
&\partial_{t}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}+(v_{1,j_{1}}ik_{x}+v_{2,j_{2}}ik_{y})\hat{f}_{k_{x},k_{y},j_{1},j_{2}}+\widehat{({\cal
F}{\cal D}f)}_{k_{x},k_{y},j_{1},j_{2}}=0,\\\
&\partial_{t}\hat{B}_{k_{x},k_{y}}=ik_{y}\hat{E}_{1,k_{x},k_{y}}-ik_{x}\hat{E}_{2,k_{x},k_{y}},\\\
&\partial_{t}\hat{E}_{1,k_{x},k_{y}}=ik_{y}\hat{B}_{k_{x},k_{y}}-\sum_{j_{1},j_{2}}v_{1,j_{1}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}+\bar{J}_{1},\\\
&\partial_{t}\hat{E}_{2,k_{x},k_{y}}=-ik_{x}\hat{B}_{k_{x},k_{y}}-\sum_{j_{1},j_{2}}v_{2,j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}+\bar{J}_{2},\end{aligned}\right.$ (72)

with the initial conditions
$\hat{f}_{k_{x},k_{y},j_{1},j_{2}}(0),\hat{B}_{k_{x},k_{y}}(0),\hat{E}_{1,k_{x},k_{y}}(0),\hat{E}_{2,k_{x},k_{y}}(0)$
satisfying the Poisson equation
$ik_{x}\hat{E}_{1,k_{x},k_{y}}(0)+ik_{y}\hat{E}_{2,k_{x},k_{y}}(0)=\sum_{j_{1},j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}(0)\Delta
v_{1}\Delta v_{2}$ for $(k_{x},k_{y})\neq(0,0)$.

For the semi-discretized system (72), we have a similar proposition as Prop
(6) in this Fourier case.

######  Proposition 7.

The semi-discretized system (72) satisfied by
$(\hat{f}_{k_{x},k_{y},j_{1},j_{2}},\hat{B}_{k_{x},k_{y}},\hat{E}_{1,k_{x},k_{y}},\hat{E}_{2,k_{x},k_{y}})(t)$
preserves the following discretized Poisson equation

$ik_{x}\hat{E}_{1,k_{x},k_{y}}(t)+ik_{y}\hat{E}_{2,k_{x},k_{y}}(t)=\sum_{j_{1},j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}(t)\Delta
v_{1}\Delta v_{2},\;\;\;(k_{x},k_{y})\neq(0,0),\forall t>0,$

provided it is satisfied at time $t=0$.

###### Proof.

As in the proof of Prop (6), we take the derivative with respect to time of
the left hand side to get

$\displaystyle\partial_{t}(ik_{x}\hat{E}_{1,k_{x},k_{y}}+ik_{y}\hat{E}_{2,k_{x},k_{y}})$
$\displaystyle=$ $\displaystyle
ik_{x}\Big{(}ik_{y}\hat{B}_{k_{x},k_{y}}-\sum_{j_{1},j_{2}}v_{1,j_{1}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}\Big{)}$
$\displaystyle+ik_{y}\Big{(}-ik_{x}\hat{B}_{k_{x},k_{y}}-\sum_{j_{1},j_{2}}v_{2,j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}\Big{)}$ $\displaystyle=$
$\displaystyle\sum_{j_{1},j_{2}}\Delta v_{1}\Delta
v_{2}\Big{(}\partial_{t}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}+\widehat{({\cal
F}{\cal D}f)}_{k_{x},k_{y},j_{1},j_{2}}\Big{)}$ $\displaystyle=$
$\displaystyle\partial_{t}\Big{(}\sum_{j_{1},j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}\Big{)},$

where we used the summation on $j_{1},j_{2}$ of the Vlasov equation together
with the fact that ${\cal D}$ is a conservative finite difference operator.
Finally, integrating the result in time between $0$ and $t$ and assuming the
relation holds at $t=0$ ends the proof. ∎

To derive a fully discrete scheme, we introduce as previously the vector
${\hat{\bf f}}_{k_{x},k_{y}}\in\mathbb{C}^{N_{v_{1}}N_{v_{2}}}$ and denote
$U(t)=({\hat{\bf
f}},{\hat{B}},{\hat{E}_{2}},{\hat{E}_{1}})_{k_{x},k_{y}}(t)\in\mathbb{M}_{N_{v_{1}}N_{v_{2}}+3,N_{v_{1}}N_{v_{2}}+3}(\mathbb{C})$,
then the system (72) can be rewritten as

$\partial_{t}U=LU+N(U),$ (73)

with

$L=\left(\begin{array}[]{lllllllccc}-i(\\!k_{x}\mbox{diag}(\vec{v_{1}})\\!\\!+\\!\\!k_{y}v_{2,1}\\!)&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\ldots&{\bf 0}_{N_{v_{1}},N_{v_{1}}}&{\bf
0}_{N_{v_{1}},1}&\dots&{\bf 0}_{N_{v_{1}},1}\\\ {\bf
0}_{N_{v_{1}},N_{v_{1}}}&-i(\\!k_{x}\mbox{diag}(\vec{v_{1}})\\!\\!+\\!\\!k_{y}v_{2,2}\\!)&\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\vdots&\dots&\vdots\\\ \vdots&\ddots&\ddots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\vdots&\dots&\vdots\\\ {\bf
0}_{N_{v_{1}},N_{v_{1}}}&\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&-i(\\!k_{x}\mbox{diag}(\vec{v_{1}})\\!\\!+\\!\\!k_{y}v_{2,N_{v_{2}}}\\!)&{\bf
0}_{N_{v_{1}},1}&\dots&{\bf 0}_{N_{v_{1}},1}\\\ {\bf
0}_{1,N_{v_{1}}}&\ldots&{\bf 0}_{1,N_{v_{1}}}&{\bf
0}_{1,N_{v_{1}}}&0&-ik_{x}&ik_{y}\\\ -\Delta v_{1}\Delta v_{2}v_{2,1}{\bf
1}&-\Delta v_{1}\Delta v_{2}v_{2,2}{\bf 1}&\ldots&-\Delta v_{1}\Delta
v_{2}v_{2,N_{v_{y}}}{\bf 1}&-ik_{x}&0&0\\\ -\Delta v_{1}\Delta
v_{2}\vec{v_{1}}&-\Delta v_{1}\Delta v_{2}\vec{v_{1}}&\ldots&-\Delta
v_{1}\Delta v_{2}\vec{v_{1}}&ik_{y}&0&0\end{array}\right)$ (74)

and

$\displaystyle\hat{\bf
f}_{k_{x},k_{y}}=\left(\begin{array}[]{lllllll}\hat{f}_{k_{x},k_{y},1,1}\\\
\vdots\\\ \hat{f}_{k_{x},k_{y},N_{v_{1}},1}\\\ \vdots\\\
\hat{f}_{k_{x},k_{y},N_{v_{1}},N_{v_{2}}}\end{array}\right)\in\mathbb{M}_{N_{v_{1}}\\!N_{v_{2}},N_{v_{1}}\\!N_{v_{2}}}(\mathbb{C})\;\;\mbox{
and }\;\;N(U)=$ $\displaystyle\left(\begin{array}[]{lllllll}(\widehat{{\cal
F{\cal D}{\bf f}}})_{k_{x},k_{y}}\\\ 0\\\ 0\\\
0\end{array}\right)\in\mathbb{M}_{N_{v_{1}}\\!N_{v_{2}}\\!+3,N_{v_{1}}\\!N_{v_{2}}\\!+3}(\mathbb{C}),$

where we denote $\vec{v_{1}}$ the vector with components
$(\vec{v_{1}})_{j_{1}}=v_{1,\min}+j_{1}\Delta v_{1}$ and
${v_{2,j_{2}}}=v_{2,\min}+j_{2}\Delta v_{2}$. In the same spirit as
previously, diag$(\vec{v_{1}})$ denotes the diagonal matrix whose diagonal is
composed of $\vec{v_{1}}$.

We now study the time discretization of (73) based on a Lawson scheme which
requires to know $\exp(tL)$ with $L$ given above. Similar (but more tedious)
calculations to those performed in the proof of Prop 4 enable to get an
explicit expression of $\exp(tL)$. To end this section, we prove that the
following scheme, with $U^{n}\approx U(t^{n}),t^{n}=n\Delta t,\Delta t>0$ and
the notations introduced above

$U^{n+1}=\exp(\Delta tL)U^{n}+\Delta t\exp(\Delta tL)N(U^{n}),$ (76)

that approximates the ODE (73) preserves a discrete Poisson equation. This is
the object of the following proposition.

######  Proposition 8.

The Lawson scheme (76) satisfied by $U^{n}=(\hat{\bf
f},\hat{B},\hat{E}_{2},\hat{E}_{1})_{k_{x},k_{y}}^{n}$ preserves the following
Poisson equation

$ik_{x}\hat{E}_{1,k_{x},k_{y}}^{n}+ik_{y}\hat{E}_{2,k_{x},k_{y}}^{n}=\sum_{j_{1},j_{2}}\hat{f}^{n}_{k_{x},k_{y},j_{1},j_{2}}\Delta
v_{1}\Delta v_{2},\;\;(k_{x},k_{y})\neq(0,0),\forall n\in\mathbb{N}^{\star},$

provided that it is satisfied at the initial time $n=0$.

###### Proof.

First, we need to know the shape of $\exp(\Delta tL)$. From the one-
dimensional calculations and from [4], we have

$\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!e^{L\Delta
t}\\!\\!=\\!\\!\\!\left(\begin{array}[]{lllllllccc}e^{-i\Delta
t(k_{x}\vec{v_{1}}+k_{y}v_{2,1})}&\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}\ldots&\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!{\bf 0}_{N_{v_{1}},1}&\\!\\!\\!{\bf
0}_{N_{v_{1}},1}&\\!\\!\\!{\bf 0}_{N_{v_{1}},1}\\\ {\bf
0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!\\!\\!\\!e^{-i\Delta
t(k_{x}\vec{v_{1}}+k_{y}v_{2,2})}&\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}\ldots&\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!\vdots&\\!\\!\\!\vdots&\\!\\!\\!\vdots\\\
\vdots&\\!\\!\\!\\!\\!\\!\ddots&\\!\\!\\!\ddots&\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!\vdots&\\!\\!\\!\vdots&\\!\\!\\!\vdots\\\
{\bf 0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!\\!\\!\\!\ldots&\\!\\!\\!\\!\\!\\!{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\\!\\!\\!\\!\\!\\!e^{-i\Delta
t(k_{x}\vec{v_{1}}+k_{y}v_{2,N_{v_{2}}})}&\\!\\!\\!{\bf 0}_{N_{v_{1}},1}&{\bf
0}_{N_{v_{1}},1}&{\bf 0}_{N_{v_{1}},1}\\\ {}_{e}{\cal
B}_{k_{x},k_{y},\star,1}&\\!\\!\\!\\!\\!\\!{}_{e}{\cal
B}_{k_{x},k_{y},\star,2}&\\!\\!\\!\\!\\!\\!\ldots&\\!\\!\\!\\!\\!\\!{}_{e}{\cal
B}_{k_{x},k_{y},\star,N_{v_{2}}}&\\!\\!\\!\\!\\!\\!\cos(|k|\Delta
t)&\\!\\!\\!\\!\\!\\!\frac{-ik_{x}\\!\sin(|k|\Delta
t)}{|k|}&\\!\\!\\!\\!\\!\\!\frac{ik_{y}\\!\sin(|k|\Delta
t)}{|k|}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\\ {}_{e}{\cal
E}_{2,k_{x},k_{y},\star,1}&\\!\\!\\!\\!\\!\\!{}_{e}{\cal
E}_{2,k_{x},k_{y},\star,2}&\\!\\!\\!\\!\\!\\!\ldots&\\!\\!\\!\\!\\!\\!{}_{e}{\cal
E}_{2,k_{x},k_{y},\star,N_{v_{2}}}&\\!\\!\\!\\!\\!\\!\\!\frac{-ik_{x}\\!\sin(|k|\Delta
t)}{|k|}&\\!\\!\\!\\!\\!\frac{k_{x}^{2}\\!\cos(|k|\Delta
t)+k_{y}^{2}}{|k|^{2}}&\\!\\!\\!\\!\\!\\!\frac{k_{x}k_{y}\\!(1\\!-\\!\cos(\Delta
t|k|)}{|k|^{2}}\\\ {}_{e}{\cal
E}_{1,k_{x},k_{y},\star,1}&\\!\\!\\!\\!\\!\\!{}_{e}{\cal
E}_{1,k_{x},k_{y},\star,2}&\\!\\!\\!\\!\\!\\!\ldots&\\!\\!\\!\\!\\!\\!{}_{e}{\cal
E}_{1,k_{x},k_{y},\star,N_{v_{2}}}&\\!\\!\\!\\!\\!\\!\\!\frac{ik_{y}\\!\sin(|k|\Delta
t)}{|k|}&\\!\\!\\!\\!\\!\\!\frac{k_{x}k_{y}(1-\cos(\Delta
t|k|)}{|k|^{2}}&\\!\\!\\!\\!\\!\frac{k_{y}^{2}\cos(\Delta
t|k|)+k_{x}^{2}}{|k|^{2}}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\end{array}\right)$

where we used the fact that the exponential of the homogeneous Maxwell part is

$\exp\left(t\left(\begin{array}[]{lllllllccc}0&-ik_{x}&ik_{y}\\\
-ik_{x}&0&0\\\
ik_{y}&0&0\end{array}\right)\right)=\left(\begin{array}[]{lllllllccc}\cos(t|k|)&-\frac{ik_{x}\sin(t|k|)}{|k|}&\frac{ik_{y}\sin(t|k|)}{|k|}\\\
-\frac{ik_{x}\sin(t|k|)}{|k|}&\frac{k_{x}^{2}\cos(t|k|)+k^{2}_{y}}{|k|^{2}}&\frac{k_{x}k_{y}(1-\cos(t|k|))}{|k|^{2}}\\\
\frac{ik_{y}\sin(t|k|)}{|k|}&\frac{k_{x}k_{y}(1-\cos(t|k|))}{|k|^{2}}&\frac{k_{y}^{2}\cos(t|k|)+k_{x}^{2}}{|k|^{2}}\end{array}\right),$

and the vectors $({}_{e}{\cal B},{}_{e}{\cal E}_{2},{}_{e}{\cal
E}_{1})_{k_{x},k_{y},\star,j_{2}}\in\mathbb{\cal C}^{3N_{v_{1}}}$ for all
$k_{x},k_{y}$ and $j_{2}=1,\dots,N_{v_{2}}$ will be given below. From the
components on $\hat{\bf f}$, since we get a diagonal matrix, we have

$\hat{f}_{k_{x},k_{y},j_{1},j_{2}}(t)=e^{-i(t-t^{n})(k_{x}v_{1,j_{1}}+k_{y}v_{2,j_{2}})}f_{k_{x},k_{y},j_{1},j_{2}}^{n}$
(77)

which can be inserted in the Maxwell part to compute the vectors $({}_{e}{\cal
B},{}_{e}{\cal E}_{2},{}_{e}{\cal E}_{1})_{k_{x},k_{y},\star,j_{2}}$. To do
so, we consider the vector $V(t)=(\hat{B},\hat{E}_{2},\hat{E}_{1})^{T}(t)$
which solves the following ODE

$\partial_{t}V={\cal J}V+R\mbox{ with }R(t)=\left(\begin{array}[]{lll}0\\\
\sum_{j_{1},j_{2}}v_{2,j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}(t)\Delta
v_{1}\Delta v_{2}\\\
\sum_{j_{1},j_{2}}v_{1,j_{1}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}(t)\Delta
v_{1}\Delta v_{2}\end{array}\right)\mbox{ and }{\cal
J}=\left(\begin{array}[]{lllllllccc}0&-ik_{x}&ik_{y}\\\ -ik_{x}&0&0\\\
ik_{y}&0&0\end{array}\right).$

The variation of constant formula gives

$V(t^{n+1})=e^{\Delta t{\cal
J}}V(t^{n})+\int_{t^{n}}^{t^{n+1}}e^{(t^{n+1}-t){\cal J}}R(t)dt.$

The calculations for the integral term involve the following integral term

$\displaystyle{\cal I}_{1}(\hat{\bf f}^{n})$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\Big{(}\frac{-ik_{x}\sin(|k|(t^{n+1}-t))}{|k|}R_{2}(t)+\frac{ik_{y}\sin(|k|(t^{n+1}-t)))}{|k|}R_{3}(t)\Big{)}dt$
$\displaystyle=$
$\displaystyle\sum_{j_{1},j_{2}}\Big{(}\frac{(-ik_{x}v_{2,j_{2}}+ik_{y}v_{1,j_{1}})}{|k|}{\cal
C}\Big{)}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}^{n}\Delta v_{1}\Delta v_{2},$
$\displaystyle=$ $\displaystyle\sum_{j_{1},j_{2}}{}_{e}{\cal
B}_{k_{x},k_{y},j_{1},j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}^{n},$
$\displaystyle{\cal I}_{2}(\hat{\bf f}^{n})$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\Big{(}\frac{k^{2}_{x}\cos(|k|(t^{n+1}-t))+k_{y}^{2}}{|k|^{2}}R_{2}(t)+\frac{k_{x}k_{y}(1-\cos(|k|(t^{n+1}-t)))}{|k|^{2}}R_{3}(t)\Big{)}dt$
$\displaystyle=$
$\displaystyle\sum_{j_{1},j_{2}}\Big{(}v_{2,j_{2}}\frac{k^{2}_{x}{\cal
A}+k_{y}^{2}{\cal B}}{|k|^{2}}+v_{1,j_{1}}\frac{k_{x}k_{y}({\cal B}-{\cal
A})}{|k|^{2}}\Big{)}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}^{n}\Delta v_{1}\Delta
v_{2},$ $\displaystyle=$ $\displaystyle\sum_{j_{1},j_{2}}{}_{e}{\cal
E}_{2,k_{x},k_{y},j_{1},j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}^{n},$
$\displaystyle{\cal I}_{3}(\hat{\bf f}^{n})$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\Big{(}\frac{k_{x}k_{y}(1-\cos(|k|(t^{n+1}-t)))}{|k|^{2}}R_{2}(t)+\frac{k^{2}_{y}\cos(|k|(t^{n+1}-t))+k_{x}^{2}}{|k|^{2}}R_{3}(t)\Big{)}dt$
$\displaystyle=$
$\displaystyle\sum_{j_{1},j_{2}}\Big{(}v_{2,j_{2}}\frac{k_{x}k_{y}({\cal
B}-{\cal A})}{|k|^{2}}+v_{1,j_{1}}\frac{k_{y}^{2}{\cal A}+k^{2}_{x}{\cal
B}}{|k|^{2}}\Big{)}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}^{n}\Delta v_{1}\Delta
v_{2},$ $\displaystyle=$ $\displaystyle\sum_{j_{1},j_{2}}{}_{e}{\cal
E}_{1,k_{x},k_{y},j_{1},j_{2}}\hat{f}_{k_{x},k_{y},j_{1},j_{2}}^{n},$

where the time integrals are

$\displaystyle{\cal A}$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\cos(|k|(t^{n+1}-t))e^{-i(t-t^{n})(k\cdot
v_{j})}dt,$ $\displaystyle{\cal B}$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}e^{-i(t-t^{n})(k\cdot
v_{j})}dt=\frac{1}{ik\cdot v_{j}}(1-e^{-i\Delta tk\cdot v_{j}}),$
$\displaystyle{\cal C}$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\sin(|k|(t^{n+1}-t))e^{-i(t-t^{n})(k\cdot
v_{j})}dt,$

with $k\cdot v_{j}=k_{x}v_{1,j_{1}}+k_{y}v_{2,j_{2}}$. To check the
conservation of the Poisson equation, one focuses on the equations on
$\hat{E}_{2}$ and $\hat{E}_{1}$ only. Thanks to the above calculations, we can
write down the update of $\hat{E}_{2},\hat{E}_{1}$ using the first order
Lawson scheme

$\displaystyle\hat{E}_{2}^{n+1}$ $\displaystyle=$ $\displaystyle(e^{\Delta
t{\cal J}}V^{n})_{2}+{\cal I}_{2}(\hat{\bf f}^{n}-\Delta t\widehat{({\cal
F}{\cal D}{\bf f}^{n})}),$ $\displaystyle\hat{E}_{1}^{n+1}$ $\displaystyle=$
$\displaystyle(e^{\Delta t{\cal J}}V^{n})_{3}+{\cal I}_{3}(\hat{\bf
f}^{n}-\Delta t\widehat{({\cal F}{\cal D}{\bf f}^{n})}).$

Thus, it remains to compute $ik_{x}\hat{E}_{1}^{n+1}+ik_{y}\hat{E}_{2}^{n+1}$
using the last relations. First, one can check easily that

$\displaystyle ik_{y}(e^{\Delta t{\cal J}}V^{n})_{2}+ik_{x}(e^{\Delta t{\cal
J}}V^{n})_{3}$ $\displaystyle=$ $\displaystyle
ik_{y}\Big{(}\frac{-ik_{x}\sin(|k|\Delta
t)}{|k|}\hat{B}^{n}+\frac{k^{2}_{x}\cos(|k|\Delta
t)+k_{y}^{2}}{|k|^{2}}\hat{E}_{2}^{n}+\frac{k_{x}k_{y}(1-\cos(|k|\Delta
t))}{|k|^{2}}\hat{E}_{1}^{n}\Big{)}$
$\displaystyle+ik_{x}\Big{(}\frac{ik_{y}\sin(|k|\Delta
t)}{|k|}\hat{B}^{n}+\frac{k_{x}k_{y}(1-\cos(|k|\Delta
t))}{|k|^{2}}\hat{E}_{2}^{n}+\frac{k^{2}_{y}\cos(|k|\Delta
t)+k_{x}^{2}}{|k|^{2}}\hat{E}_{1}^{n}\Big{)}$ $\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}.$

Then, we have

$\displaystyle ik_{y}\hat{E}_{2}^{n+1}+ik_{x}\hat{E}_{1}^{n+1}$
$\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}+ik_{y}{\cal
I}_{2}(\hat{f}^{n}-\Delta t\widehat{({\cal F}{\cal D}f^{n})})+ik_{x}{\cal
I}_{3}(\hat{f}^{n}-\Delta t\widehat{({\cal F}{\cal D}f^{n})})$
$\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}+ik_{y}\sum_{j_{1},j_{2}}v_{2,j_{2}}(\hat{f}^{n}-\Delta
t\widehat{({\cal F}{\cal D}f^{n})}){\cal B}\Delta v_{2}\Delta v_{1}$
$\displaystyle+ik_{x}\sum_{j_{1},j_{2}}v_{1,j_{1}}(\hat{f}^{n}-\Delta
t\widehat{({\cal F}{\cal D}f^{n})}){\cal B}\Delta v_{2}\Delta v_{1}$
$\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}+\Delta v_{2}\Delta
v_{1}\sum_{j_{1},j_{2}}(ik\cdot v_{j})(\hat{f}^{n}-\Delta t\widehat{({\cal
F}{\cal D}f^{n})})\frac{1}{ik\cdot v_{j}}(1-e^{-i\Delta tk\cdot v_{j}})$
$\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}+\Delta v_{2}\Delta
v_{1}\Big{[}\sum_{j_{1},j_{2}}(1-e^{-i\Delta tk\cdot v_{j}})f^{n}$
$\displaystyle-\Delta t\sum_{j_{1},j_{2}}\widehat{({\cal F}{\cal
D}f^{n}})+\Delta t\sum_{j_{1},j_{2}}e^{-i\Delta tk\cdot v_{j}}\widehat{({\cal
F}{\cal D}f^{n}})\Big{]}$ $\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}+\Delta v_{2}\Delta
v_{1}\Big{[}\sum_{j_{1},j_{2}}(1-e^{-i\Delta tk\cdot
v_{j}})\hat{f}^{n}+\sum_{j_{1},j_{2}}e^{-i\Delta tk\cdot
v_{j}}\hat{f}^{n}-\hat{f}^{n+1}\Big{]}$ $\displaystyle=$ $\displaystyle
ik_{y}\hat{E}_{2}^{n}+ik_{x}\hat{E}_{1}^{n}+\Delta v_{2}\Delta
v_{1}\sum_{j_{1},j_{2}}(\hat{f}^{n}-\hat{f}^{n+1}),$

where we used the update for $f$: $\hat{f}^{n+1}=e^{-i\Delta tk\cdot
v_{j}}\hat{f}^{n}-\Delta te^{-i\Delta tk\cdot v_{j}}\widehat{({\cal F}{\cal
D}f^{n}})$ and the conservation property of the discrete operator ${\cal D}$.
Then, if the Poisson equation is satisfied at iteration $n$, it is propagated
to the next iteration, which concludes the proof. ∎

##  5 Numerical experiments

In this section, we perform numerical experiments for linear transport
problems and Vlasov equations. First, we study the different order of
convergence on a linear problem and then, we present some numerical results of
the exponential DG solutions for Vlasov equations in 1dx-1dv and 1dx-2dv
cases.

###  5.1 2D linear passive-transport problems

We consider the following two-dimension linear transport equation

$\partial_{t}u+\partial_{x}u+\partial_{v}u=0,\ (x,v)\in[0,2\pi]^{2}$ (78)

with the initial condition $u(x,v,0)=\sin(x+v)$ and periodic boundary
condition. The exact solution is $u(x,v,t)=\sin(x+v-2t)$ which enables us to
check the different order of convergence. Indeed, for a Lawson scheme based on
a underlying Runge-Kutta method RK($m,s)$ (order $m$, $s$ stages, a DG space
approximation with $P_{k}$ and a finite difference approximation in $v$ of
order $4$ (which means ${\cal D}$ is chosen as a 4th order centered finite
difference operator CD4), we expect the following estimate

$\|{\bf u}(t^{n})-{\bf u}_{h}^{n}\|_{L^{2}}\leq C(\Delta x^{k}+\Delta
v^{4}+\Delta t^{m})$

Here we use the 3rd order Lawson-RK method for the time discretization, with a
final time $T=1$, and consider different parameters to test the convergence
rates in $t,x$ and $v$. We firstly take $N_{v}=320,\Delta t=0.01$, and
consider different mesh size $N_{x}$ to check the convergence rate of DG in x
direction for both central and upwind fluxes (in this linear case, upwind
fluxes can be considereed easily). Table 1 shows the $L^{\infty}$ and $L^{2}$
errors, the associated orders of convergence for DG-$P^{k}$ for $k=1,2$ in
$x-$ direction. The optimal convergence rate for DG is clearly obtained. In
particular, the sub-optimal and optimal rates are observed according to the
choice of the flux and to the oddness of $k$, as discussed in [23]. Then, we
study the convergence in $v-$ direction. We take $k=5,N_{x}=32,\Delta t=0.01$,
and consider different mesh size $N_{v}$ to check the convergence rate of the
fourth order approximation of ${\cal D}$ (CD4). Table 2 shows the expected
convergence (note that only central fluxes are considered in this case).
Finally, to check time accuracy, we take $k=5,N_{x}=16,N_{v}=32$ and $\Delta
t=0.0001$ to compute a reference solution. Then we get the error table for
different time step sizes $\Delta t$. From Table 3, the expected 3rd order
convergence is observed for all cases.

| $N_{x}$ | central flux | upwind flux
---|---|---|---
$L^{\infty}$-error | order |  $L^{2}$-error | order |  $L^{\infty}$-error | order |  $L^{2}$-error | order
$P^{1}$ | $10$ | 1.62E-01 | - | 3.76E-01 | - | 4.78E-02 | - | 9.51E-02 | -
$20$ | 7.66E-02 | 1.08 | 1.86E-01 | 1.02 | 1.27E-02 | 1.91 | 2.42E-02 | 1.98
$40$ | 3.70E-02 | 1.05 | 9.24E-02 | 1.01 | 3.25E-03 | 1.97 | 6.07E-03 | 1.99
$80$ | 1.82E-02 | 1.03 | 4.60E-02 | 1.00 | 8.23E-04 | 1.99 | 1.52E-03 | 2.00
| $160$ | 9.01E-03 | 1.01 | 2.30E-02 | 1.00 | 2.07E-04 | 1.99 | 3.79E-04 | 2.00
$P^{2}$ | $10$ | 2.29E-03 | - | 3.76E-03 | - | 2.52E-03 | - | 4.67E-03 | -
$20$ | 2.65E-04 | 3.11 | 4.53E-04 | 3.05 | 3.07E-04 | 3.04 | 5.83E-04 | 3.00
$40$ | 3.25E-05 | 3.02 | 5.63E-05 | 3.01 | 3.84E-05 | 3.00 | 7.29E-05 | 3.00
$80$ | 4.05E-06 | 3.01 | 7.03E-06 | 3.00 | 4.79E-06 | 3.00 | 9.11E-06 | 3.00
| $160$ | 5.12E-07 | 2.98 | 8.99E-07 | 2.97 | 5.95E-07 | 3.00 | 1.16E-06 | 2.98
Table 1: Linear transport equation: $L^{\infty}$ and $L^{2}$-norm space errors of the Lawson-DG scheme with $P^{1}$, $P^{2}$ (CD4 in velocity ($N_{v}=320$) and RK(3,3) in time ($\Delta t=0.01$)).  | $N_{v}$ | central flux
---|---|---
$L^{\infty}$-error | order |  $L^{2}$-error | order
CD4 | $8$ | 1.18E-02 | - | 5.24E-02 | -
$16$ | 7.78E-04 | 3.92 | 3.50E-03 | 3.92
$32$ | 4.93E-05 | 3.98 | 2.19E-04 | 3.98
$64$ | 3.09E-06 | 4.00 | 1.37-05 | 4.00
| $128$ | 1.98E-07 | 3.97 | 8.78E-07 | 3.97
Table 2: Linear transport equation: $L^{\infty}$ and $L^{2}$-norm velocity errors of the Lawson-DG scheme with CD4 (DG-$P^{5}$ in space ($N_{x}=32$) and RK(3,3) ($\Delta t=0.01$)). | $\Delta t$ | central flux
---|---|---
$L^{\infty}$-error | order |  $L^{2}$-error | order
RK(3,3) | $0.1000$ | 4.17E-05 | - | 1.85E-04 | -
$0.0500$ | 5.21E-06 | 2.99 | 2.31E-05 | 2.99
$0.0250$ | 6.51E-07 | 3.00 | 2.89E-06 | 3.00
$0.0125$ | 8.14E-08 | 3.00 | 3.61E-07 | 3.00
| $0.00625$ | 1.02E-08 | 3.00 | 4.52E-08 | 3.00
Table 3: Linear transport equation: $L^{\infty}$ and $L^{2}$-norm time errors
of the Lawson-DG scheme with RK(3,3) (DG-$P^{5}$ in space ($N_{x}=16$) and CD4
($N_{v}=32$)).

###  5.2 Vlasov-Ampère equation

We firstly consider the following initial condition for Landau damping

$f_{0}(x,v)=\frac{1}{\sqrt{2\pi}}e^{-v^{2}/2}(1+\alpha\cos(kx)),$

where we take $x\in[0,2\pi/k],\ k=0.5,\ v\in[-9,9]$ and $\alpha=10^{-3}$. Here
we still use DG method for space discretization, the finite difference method
$CD4$ in $v$ direction and the 3rd Lawson-RK method for time discretization
(see Section 3.1). The numerical parameters are chosen as follows: $\Delta
t=0.1$, $N_{x}=31(P^{2}),N_{v}=121$.

In Figure 1, the time evolution of the electric energy $\|E(t)\|_{L^{2}}$ is
displayed in semi-log scale (with the corresponding damping rate in red) and
the deviation of the total energy ${\cal H}(t)-{\cal H}(0)$ with ${\cal
H}(t)=\int v^{2}f(t,x,v)dxdv+\int E^{2}(t,x)dx$. The expected behaviors
(correct damping rate and good energy conservation) are recovered.

Figure 1: Vlasov-Ampère equation (Landau damping): time evolution of the
electric energy in semi-log scale (left) and of the deviation of the total
energy (right). Lawson-DG RK(3,3) and $P^{2}$ ($\Delta t=0.1$,
$N_{x}=31(P^{2}),N_{v}=121$).

We consider a second test called the two stream instability test with the
initial condition

$f_{0}(x,v)=\frac{1}{\sqrt{2\pi}}v^{2}e^{-v^{2}/2}(1+\alpha\cos(kx)),$

for which the same physical and numerical parameters as previously are kept
except the final time which is $T=300$. In Figures 2, we plot the time
evolution of the electric energy in semi-log scale (and the corresponding
instability rate in red) and the deviation of the total energy. For this test,
a linear instability is first observed (up to $t\approx 30$) during which a
vortex in phase space is created (see 2), and it is followed by a nonlinear
phase. These two behaviors are well reproduced by the scheme even if the mesh
is quite coarse. In partiular, even if the vortex is well captured, we can
observe spurious oscillations due to the use of central schemes. Note that the
Poisson equation is satisfied in both cases up to machine accuracy.

Figure 2: Vlasov-Ampère equation (two stream instability): time evolution of
the electric energy in semi-log scale (top left) and of the deviation of the
total energy (top right), snapshot of $f(t=0)$ (bottom left) and snapshot of
$f(t=300)$ (bottom right). Lawson-DG RK(3,3) and $P^{2}$ ($\Delta t=0.1$,
$N_{x}=31(P^{2}),N_{v}=121$).

###  5.3 Vlasov-Maxwell equations 1dx-2dv

We consider the Weibel instability [32] by consideering the Vlasov-Maxwell
1dx-2dv model studied in Section 3.2 with the initial distribution and fields
are of the form

$f(t=0,x,v_{1},v_{2})=\frac{1}{2\pi\sigma_{1}\sigma_{2}}\exp\left(-\frac{1}{2}\left(\frac{v_{1}^{2}}{\sigma_{1}^{2}}+\frac{v_{2}^{2}}{\sigma_{2}^{2}}\right)\right)(1+\alpha\cos(k_{w}x)),\
x\in[0,2\pi/k),$ $B(t=0x)=\beta\cos(kx),\;\;\;E_{2}(t=0,x)=0,$

and $E_{1}(x,t=0)$ is imposed from the Poisson equation. We choose the
parameters $\sigma_{1}=0.02/\sqrt{2},\ \sigma_{2}=\sqrt{12}\sigma_{1},\
k=1.25,\ \alpha=0,\ \beta=-10^{-4}$ for our test, which gives a growth rate of
0.02784 by solving the dispersion relation (see Weibel [32]). For the
numerical simulations up to a final time $T=500$, we still use DG method for
space discretization in $x-$ direction, finite difference method with $CD4$ in
$v$ direction and Lawson-RK(3,3) method for time discretization (Lawson-
RK(3,3)-DG CD4) and consider $\Delta t=0.5$,
$N_{x}=21(P^{2}),N_{v_{1}}=N_{v_{2}}=44$. For comparison, we also consider
Fourier method for space discretization in $x-$ direction, finite difference
with a third order upwind (UP3) in $v$ direction and Lawson-RK(3,3) method in
time (Lawson-RK(3,3)-Fourier UP3) with $\Delta t=0.5$,
$N_{x}=64,N_{v_{1}}=N_{v_{2}}=88$.

In Figure 3, we show the time evolution of the electromagnetic energies
$\|B(t)\|_{L^{2}}$, $\|E_{y}(t)\|_{L^{2}}$, $\|E_{x}(t)\|_{L^{2}}$ (in semi-
log scale) obtained by the two methods. First, we can observe that the
theoretical growth rate is in very good agreement with the two numerical
solution. Second, the two methods are very close up to time $t\approx 300$
(which corresponds to the end of the linear phase) and slightly differs for
larger times. We also show the evolution of the relative total energy $|{\cal
H}(t)-{\cal H}(0)|/{\cal H}(0)$ in Figure 4 for Lawson-RK(3,3)-DG CD4 scheme
without (left) and with (right) the energy correction step presented in [4] in
the Vlasov-Ampère case. This projection approach enables to modify the unknown
by a suitable coefficient which is of order the scheme so that the total
energy is preserved almost up to machine error $10^{-12}$.

Figure 3: Vlasov-Maxwell equation (Weibel instability): time evolution of the
electromagnetic energies in semi-log scale together with the analytic growth
rate. Left: Lawson-RK(3,3)-DG CD4. Right: Lawson-RK(3,3)-Fourier UP3.

Figure 4: Vlasov-Maxwell equation (Weibel instability): time evolution the
relative total energy of Lawson-RK(3,3)-DG CD4 without (left) and with (right)
energy correction step.

The second test for Vlasov-Maxwell equation we considered is the streaming
Weibel instability [7, 9] for which the initial condition is

$f(t=0,x,v_{1},v_{2})=\frac{1}{2\pi\sigma^{2}}\exp\left(-\frac{v_{1}^{2}}{2\sigma^{2}}\right)\left(\delta\exp\left(-\frac{(v_{2}-v_{0,1})^{2}}{2\sigma^{2}}\right)+(1-\delta)\exp\left(-\frac{(v_{2}-v_{0,2})^{2}}{2\sigma^{2}}\right)\right),$
$B(t=0,x,t=0)=\beta\cos(kx),\;\;E_{2}(t=0,x)=0,$

and $E_{1}(t=0,x)=0$ from the Poisson equation. We choose the parameters
$\sigma=0.1/\sqrt{2},\ k=0.2,\ \beta=-10^{-3},\ v_{0,1}=0.5,\ v_{0,2}=-0.1$
and $\delta=1/6$ for our test. The growth rate of $E_{2}$ is 0.03 [7]. For the
two schemes Lawson-RK(3,3)-DG CD4 and Lawson-RK(3,3)-Fourier UP3, we take the
same parameters as in the previous case but $\Delta t=0.1$ for stability
reasons. We show the result in Figure 5 in which the time evolution of the
$L^{2}$ norm of the electromagnetic fields are displayed. First, we observe a
good agreement with the theoretical growth rate for these two schemes and some
deviations in the nonlinear phase can be observed. For this case, we also
consider the correction on the total energy and plot the time history of the
relative total energy in Figure 6 for Lawson-RK(3,3)-DG CD4 without (left) and
with (right) energy correction step. The total energy can be well preserved
for exponential DG approximations if it is stable, and the relative error is
greatly improved with energy correction step without affecting the accuracy of
the scheme. Note that the Poisson equation is satisfied in both cases up to
machine accuracy.

Figure 5: Vlasov-Maxwell equation (streaming Weibel instability): time
evolution of the electromagnetic energies in semi-log scale together with the
analytic growth rate. Left: Lawson-RK(3,3)-DG CD4. Right: Lawson-
RK(3,3)-Fourier UP3.

Figure 6: Vlasov-Maxwell equation (streaming Weibel instability): time
evolution the relative total energy of Lawson-RK(3,3)-DG CD4 without (left)
and with (right) energy correction step.

##  6 Conclusion

In this work, we constructed and implemented a new exponential DG method for
Vlasov equation, extending the previous works [14, 12, 4] on this topic where
Fourier method in space were used. These methods allow to derive high order
accuracy in time, space and velocity, still ensuring stability without the
restrictive CFL type constraint coming from the linear part. Moreover, a
discrete Poisson equation is satisfied and a projection technique enables to
preserve the total energy. The extension to DG turns out to be an interesting
alternative compared to previous approach based on Fourier which is restricted
to cartesian domains with periodic boundary conditions. The approach only
involves the calculation of exponential of DG-matrices of size $(k+1)N$ with
$k$ the DG degree and $N$ the number of points in the space direction and we
observe that thanks to the specific structure, this is also true in multi-
dimensions.

One interesting extension is to perform an efficient implementation of the
method in two-dimension in space to exploit the Kronecker structure. We also
plan to investigate the extension of this approach to problems involving
boundary conditions, for which monotone fluxes are more appropriate but
requires to study the stability of the scheme.

##  Appendix A Appendix A: The DG matrix construction

The goal of this appendix is to give some details on the construction of the
DG-matrix (8) which is at the core of our scheme. To do so, we consider the
1dx transport equation (1) with coefficient $a=1$, for which we wrote in (2)
the semi-discrete DG scheme. To get the DG matrix, we consider
$\psi=\xi^{\ell}_{j}$ and
$u_{h}(t,x)=\sum_{m=0}^{k}u^{m}_{j}(t)\xi^{m}_{j}(x)$ in (2) to get

$\sum_{m=0}^{k}\Big{[}(\partial_{t}{u}_{j}^{m}(t)\xi^{m},\xi^{\ell})_{I_{j}}-({u}_{j}^{m}(t)\xi^{m},\partial_{x}\xi^{\ell})_{I_{j}}\Big{]}+[\\{u_{h}\\}\xi^{\ell}]^{x_{j+1/2}}_{x_{j-1/2}}=0,$
(79)

where $\ell=0,\dots,k,j=1,\dots,N$ and where the central flux
$\\{u_{h}\\}|_{x_{j\pm 1/2}}=(u_{h}^{+}(x_{j\pm 1/2})+u_{h}^{-}(x_{j\pm
1/2}))/2$ is considered. Thus we obtain for the boundary term

$\displaystyle[\\{u_{h}\\}\xi^{\ell}]^{x_{j+1/2}}_{x_{j-1/2}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{m=0}^{k}\left[\left({u}_{j}^{m}(t)\xi^{m,-}(x_{j+1/2})+{u}_{j+1}^{m}(t)\xi^{m,+}(x_{j+1/2})\right)\xi^{\ell,-}(x_{j+1/2})\right.$
$\displaystyle\left.-\left({u}_{j-1}^{m}(t)\xi^{m,-}(x_{j-1/2})+{u}_{j}^{m}(t)\xi^{m,+}(x_{j-1/2})\right)\xi^{l,+}(x_{j-1/2})\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{m=0}^{k}\left[\left({u}_{j}^{m}(t)(1/2)^{m}+{u}_{j+1}^{m}(t)(-1/2)^{m}\right)(1/2)^{\ell}\right.$
$\displaystyle\left.-\left({u}_{j-1}^{m}(t)(1/2)^{m}+{u}_{j}^{m}(t)(-1/2)^{m}\right)(-1/2)^{\ell}\right].$

We denote
$\mathbf{u}_{j}(t)=({u}_{j}^{0}(t),{u}_{j}^{1}(t),....,{u}_{j}^{k}(t))^{T},\
j=1,2,...,N,$ then we can rewrite the DG discretization as

$\displaystyle\Delta x\left(\begin{array}[]{llll}M&&&\\\ &M&&\\\ &&\ddots&\\\
&&&M\\\
\end{array}\right)\partial_{t}\left(\begin{array}[]{llll}\mathbf{u}_{1}\\\
\mathbf{u}_{2}\\\ \vdots\\\ \mathbf{u}_{N}\end{array}\right)$
$\displaystyle-\left(\begin{array}[]{llll}D_{1}&D_{2}&\ldots&D_{3}\\\
D_{3}&D_{1}&D_{2}&\\\ &\ddots&\ddots&D_{2}\\\ D_{2}&\ldots&D_{3}&D_{1}\\\
\end{array}\right)\left(\begin{array}[]{llll}\mathbf{u}_{1}\\\
\mathbf{u}_{2}\\\ \vdots\\\ \mathbf{u}_{N}\end{array}\right)={\bf 0},$ (80)

where $M,D_{i}\in\mathbb{M}_{(k+1),(k+1)}(\mathbb{R}),i=1,2,3$ are given by

$\displaystyle
M_{\ell,m}=\frac{(1/2)^{m+\ell-1}}{m+\ell-1}[1-(-1)^{m+\ell-1}],\\!\\!\\!\\!\\!\\!$
$\displaystyle(D_{2})_{\ell,m}=(-1)^{m}(1/2)^{m+\ell-1},$
$\displaystyle(D_{1})_{\ell,m}=(1/2)^{m+\ell-2}\Big{(}\frac{\ell-1}{m+\ell-2}-\frac{1}{2}\Big{)}[1-(-1)^{m+\ell-2}],\\!\\!\\!\\!\\!\\!$
$\displaystyle(D_{3})_{\ell,m}=(-1)^{\ell-1}(1/2)^{m+\ell-1},$

with $(D_{1})_{1,1}=0$ by convention. Then we have

$\Delta x\mathbfcal{M}\partial_{t}\mathbf{u}-\mathbfcal{D}\mathbf{u}=0,$

where $\mathbfcal{M}$ is a block diagonal mass matrix of size $N(k+1)$,
$\mathbfcal{D}$ is a block tridiagonal matrix of size $N(k+1)$, and
$\mathbf{u}$ is the vector containing the degree of freedom

$\mathbf{u}=(u_{1}^{0},u_{1}^{1},\dots,u_{1}^{k},u_{2}^{0},u_{2}^{1},\dots,u_{2}^{k},\dots,u_{N}^{0},u_{N}^{1},\dots,u_{N}^{k})^{T}\in\mathbb{R}^{(k+1)N}.$

Now we can rewrite the DG scheme as

$\partial_{t}\mathbf{u}(t)=A\mathbf{u}(t),\ \ A=(\Delta
x\mathbfcal{M})^{-1}\mathbfcal{D},$ (81)

where $A$ is a block circulant matrix

$A=\frac{1}{\Delta
x}\left(\begin{array}[]{lllll}M^{-1}D_{1}&M^{-1}D_{2}&\bf{0}&\ldots&M^{-1}D_{3}\\\
M^{-1}D_{3}&M^{-1}D_{1}&M^{-1}D_{2}&\bf{0}&\ldots\\\
\bf{0}&\ddots&\ddots&\ddots&M^{-1}D_{2}\\\
M^{-1}D_{2}&\bf{0}&\ldots&M^{-1}D_{3}&M^{-1}D_{1}\end{array}\right)=\frac{1}{\Delta
x}A_{1},$

or with the circblock notation
$A_{1}=\mbox{circblock}(M^{-1}D_{1},M^{-1}D_{2},0,0,...,M^{-1}D_{3})$. Since
$\mathbfcal{D}$ and $\mathbfcal{M}$ are independent of $\Delta x$, so is
$A_{1}$ defined by $A=\Delta x^{-1}A_{1}$.

##  Appendix B Appendix B: Stability and error estimate for semi-discrete DG
scheme.

In this appendix, we give some error estimate of the exponential-DG scheme for
the one dimensional linear advection equation (1). To do so, we first define
some notations about norms which will be used. For a given function $x\mapsto
v(x)$, we denote $\|v\|_{j}$ and $\|v\|_{\infty,j}$ as the $L^{2}$-norm and
$L^{\infty}$-norm of $v$ on $I_{j}$ ($j=1,\dots,N$) respectively. Moreover,

$\displaystyle\|v\|=(\sum\limits_{j}\|v\|^{2}_{j})^{\frac{1}{2}},$
$\displaystyle\|v\|_{\infty}=\max\limits_{j}\|v\|_{\infty,j},$
$\displaystyle\|v\|_{\Gamma}^{2}=\sum_{j}(|v_{j-1/2}^{+}|^{2}+|v_{j-\frac{1}{2}}^{-}|^{2}),$

where we express the value of $v$ on the left and right limits of the grid
point $x_{j+\frac{1}{2}}$ with $v_{j+\frac{1}{2}}^{-}$ and
$v_{j+\frac{1}{2}}^{+}$ respectively. Define the jump and the mean of $v$ at
$x_{j-\frac{1}{2}}$ as
$[v]_{j-\frac{1}{2}}=(v_{j-\frac{1}{2}}^{+}-v_{j-\frac{1}{2}}^{-})$ and
$\\{v\\}_{j-\frac{1}{2}}=(v^{+}_{j-\frac{1}{2}}+v^{-}_{j-\frac{1}{2}})/2$.

###  B.1 Notations for projections and some properties of approximation space

The inverse properties of the finite space $V_{h}$ will be used.

######  Lemma B.1.

When the mesh is regular, $\forall v\in V_{h},\ \exists\ C>0,\ s.t.$

$h^{2}\|\partial_{x}v\|^{2}+h\|v\|_{\Gamma}^{2}\leqslant C\|v\|^{2},$ (82)

where the positive constant $C$ is independent of $h$ and $v$.

Define the $L^{2}$-projection $P_{k}$ of $u$ into $V_{h}$ as follows:

$\displaystyle(P_{k}u,v_{h})$ $\displaystyle=(u,v_{h}),\forall v_{h}\in
V_{h}\mbox{ with }v_{h}(\chi_{j}(.,t))\in P^{k}([-1,1]).$

The following lemma states the error of these projections [10].

######  Lemma B.2.

Let $P_{h}^{\perp}q=q-P_{k}q$ is the projection error. For any smooth function
$q(x)$, $\exists\ c>0,$ such that

$\|P_{h}^{\perp}q\|_{D}+h\|\partial_{x}(P_{h}^{\perp}q)\|_{D}+h^{\frac{1}{2}}\|P_{h}^{\perp}q\|_{\infty,D}\leqslant
ch^{k+1}|q|_{k+1,D},$ (83) $\|P_{h}^{\perp}q\|_{\Gamma}\leqslant
ch^{k+\frac{1}{2}}\|\partial_{x}^{k+1}q\|,$ (84)

where the positive constant $c$ is not dependent on $h$, solely depending on
$q$, and $D$ may be $\Omega$ or $I_{j}$.

Furthermore, to avoid confusion with different constants, we denote a generic
positive constant by $C$, which is independent of the numerical solution and
the mesh size for our problem. But, the constant may dependent on the exact
solution and may have a different value in each occurrence. Moreover, for
problems considered in this paper, the exact solution is assumed to be smooth
with periodic or compactly supported boundary condition. Therefore, the exact
solution is always bounded.

We state the $L^{2}$ stability and $L^{2}$-norm error estimate for the scheme
and also give their proof.

######  Theorem B.1.

For semi-discrete DG scheme (2) with central flux, we have the
$L^{2}$-stability:

$\frac{d}{dt}||u_{h}||_{L^{2}}^{2}=0.$

###### Proof.

Take the test function $\psi=u_{h}$ in the semi-discrete scheme (2), we have

$\displaystyle\frac{1}{2}\frac{d}{dt}\int_{{I}_{j}}u_{h}^{2}dx$
$\displaystyle=-a\\{u_{h}\\}|_{{x}_{j+\frac{1}{2}}}u_{h}^{-}|_{{x}_{j+\frac{1}{2}}}+a\\{u_{h}\\}|_{{x}_{j-\frac{1}{2}}}u_{h}^{+}|_{{x}_{j-\frac{1}{2}}}+\int_{{I}_{j}}au_{h}(u_{h})_{x}dx,$
(85)
$\displaystyle=-a\\{u_{h}\\}|_{{x}_{j+\frac{1}{2}}}u_{h}^{-}|_{{x}_{j+\frac{1}{2}}}+a\\{u_{h}\\}|_{{x}_{j-\frac{1}{2}}}u_{h}^{+}|_{{x}_{j-\frac{1}{2}}}+\frac{a}{2}(u_{h}^{2})^{-}|_{{x}_{j+\frac{1}{2}}}-\frac{a}{2}(u_{h}^{2})^{+}|_{{x}_{j-\frac{1}{2}}}.$
$\displaystyle=-\frac{a}{2}u_{h}^{+}u_{h}^{-}|_{{x}_{j+\frac{1}{2}}}+\frac{a}{2}u_{h}^{-}u_{h}^{+}|_{{x}_{j-\frac{1}{2}}}.$

Sum over $j$ of above equation, the $L^{2}$-stability follows. ∎

######  Theorem B.2.

Let $T>0$, $u$ be the exact solution of problem (1), which is sufficiently
smooth with bounded derivatives. Assume $u_{h}$ is the DG approximation of
semi-discrete scheme (2) with the central flux and the approximation space
$V_{h}$ is the space consisting of $k$-th piecewise polynomial. Then it holds
that

$\|u(T)-u_{h}(T)\|_{L^{2}}\leqslant Ch^{k},$ (86)

where C is a positive constant independent on $\Delta x$.

###### Proof.

Denote error as $e_{u}=u_{h}-u$. Notice that the scheme (2) is still satisfied
with $u_{h}=u$. So, we have the error equation

$\displaystyle\int_{{I}_{j}}(\partial_{t}e_{u}\psi)dx=-a\\{e_{u}\\}|_{{x}_{j+\frac{1}{2}}}\psi^{-}|_{{x}_{j+\frac{1}{2}}}+a\\{e_{u}\\}|_{{x}_{j-\frac{1}{2}}}\psi^{+}|_{{x}_{j-\frac{1}{2}}}+\int_{{I}_{j}}ae_{u}\psi_{x}dx.$
(87)

Define
$e_{u}=u_{h}-u=(u_{h}-P_{k}u)-(u-P_{k}u)=\widetilde{e_{u}}-P_{k}^{\perp}u$.
Then taking $\psi=\widetilde{e_{u}}$, we have

$\displaystyle\int_{{I}_{j}}\partial_{t}(\widetilde{e_{u}})\widetilde{e_{u}}dx$
(88)
$\displaystyle=\int_{{I}_{j}}\partial_{t}(P_{k}^{\perp}u)\widetilde{e_{u}}dx-a\\{\widetilde{e_{u}}\\}|_{{x}_{j+\frac{1}{2}}}\widetilde{e_{u}}^{-}|_{{x}_{j+\frac{1}{2}}}+a\\{\widetilde{e_{u}}\\}|_{{x}_{j-\frac{1}{2}}}\widetilde{e_{u}}^{+}|_{{x}_{j-\frac{1}{2}}}+\int_{{I}_{j}}a\widetilde{e_{u}}(\widetilde{e_{u}})_{x}dx$
$\displaystyle+a\\{P_{k}^{\perp}u\\}|_{{x}_{j+\frac{1}{2}}}\widetilde{e_{u}}^{-}|_{{x}_{j+\frac{1}{2}}}-a\\{P_{k}^{\perp}u\\}|_{{x}_{j-\frac{1}{2}}}\widetilde{e_{u}}^{+}|_{{x}_{j-\frac{1}{2}}}-\int_{{I}_{j}}aP_{k}^{\perp}u(\widetilde{e_{u}})_{x}dx.$

By the the definition of the projections and some calculations, the right
terms $RHS$ of the error equation (88) become

$\displaystyle
RHS=-\frac{a}{2}\widetilde{e_{u}}^{+}\widetilde{e_{u}}^{-}|_{{x}_{j+\frac{1}{2}}}+\frac{a}{2}\widetilde{e_{u}}^{-}\widetilde{e_{u}}^{+}|_{{x}_{j-\frac{1}{2}}}+a\\{P_{k}^{\perp}u\\}|_{{x}_{j+\frac{1}{2}}}\widetilde{e_{u}}^{-}|_{{x}_{j+\frac{1}{2}}}-a\\{P_{k}^{\perp}u\\}|_{{x}_{j-\frac{1}{2}}}\widetilde{e_{u}}^{+}|_{{x}_{j-\frac{1}{2}}}.$
(89)

Sum over $j$,

$\frac{1}{2}\frac{d}{dt}||\widetilde{e_{u}}||^{2}=\sum_{j}RHS=-\sum_{j}a\\{P_{k}^{\perp}u\\}[\widetilde{e_{u}}]|_{{x}_{j+\frac{1}{2}}}.$

Furthermore, from Holder’s inequality, Lemma (B.1) and Lemma (B.2), we have

$\displaystyle\frac{1}{2}\frac{d}{dt}||\widetilde{e_{u}}||^{2}\leq\sum_{j}\left|a\\{P_{k}^{\perp}u\\}[\widetilde{e_{u}}]|_{{x}_{j+\frac{1}{2}}}\right|$
(90) $\displaystyle\leq
c||P_{k}^{\perp}u||_{\Gamma}||\widetilde{e_{u}}||_{\Gamma}\leq
ch^{k+\frac{1}{2}}h^{-\frac{1}{2}}||\widetilde{e_{u}}||$ $\displaystyle\leq
c\|\widetilde{e_{u}}\|^{2}+ch^{2k}.$

Thus by Gronwall’s inequality, the conclusion in Theorem B.2 follows. ∎

##  Appendix C Appendix C: Ker$(A)$ and projection $\Pi$

Here, the projection matrix $\Pi$ onto Ker$(A)$ with $A$ the DG-matrix (8) is
discussed. A general expression (for any $k,N$) turns out to be difficult and
we compute Ker$(A)$ (and the projection $\Pi$) for several practical cases. As
mentioned in Remark 2, there are mainly two cases, according to the oddness of
$(k+1)N$: $(i)$ if $(k+1)N$ is odd, $0$ is a simple eigenvalue of $A$ and
Ker$(A)$=Span$(u_{1})$ where $u_{1}\in\mathbb{R}^{(k+1)N}$ correspond
constants in the space $P^{k}$; $(ii)$ if $(k+1)N$ is even, $0$ is a double
eigenvalue of $A$ and Ker$(A)$=Span$(u_{1},u_{2})$ and
$u_{2}\in\mathbb{R}^{(k+1)N}$ has to be determined. This second case recalls
what happens for the second centered finite differences in which constant
vector belongs to the kernel but also the sequence $(-1)^{j}$.

  * •

odd case: We can check that

$u_{1}=\Big{[}1,\underbrace{0,\ldots,0}_{\in\mathbb{R}^{k}},1,\underbrace{0,\ldots,0}_{\in\mathbb{R}^{k}},\ldots,1,\underbrace{0,\ldots,0}_{\in\mathbb{R}^{k}}\Big{]}^{T}\in\mathbb{R}^{(k+1)N}.$

is a eigenvector of $A$ associated to the eigenvalue $0$. By $\Pi x=\langle
x,u_{1}\rangle u_{1}$, we get the expression of the matrix $\Pi$

$\Pi=[u_{1},{\bf 0},\ldots,{\bf 0},u_{1},{\bf 0},\ldots,{\bf
0},\ldots,u_{1},{\bf 0},\ldots,{\bf
0}\Big{]}\in\mathbb{M}_{(k+1)N,(k+1)N}(\mathbb{R}),$

with ${\bf 0}\in\mathbb{R}^{(k+1)N}$.

  * •

even case: in addition to $u_{1}$, we need to find a second eigenvector to
construct $\Pi$. We give below the expression of $u_{2}$ for some $k=1$ to
$k=5$

    * –

$k=0$, $u_{2}=\Big{[}0,1,0,1,\ldots,0,1\Big{]}^{T}\in\mathbb{R}^{N}.$

    * –

$k=1$,
$u_{2}=\Big{[}\underbrace{0,1}_{\in\mathbb{R}^{2}},\underbrace{0,1}_{\in\mathbb{R}^{2}},\ldots,\underbrace{0,1}_{\in\mathbb{R}^{2}}\Big{]}^{T}\in\mathbb{R}^{2N}.$

    * –

$k=2$,
$u_{2}=\Big{[}\underbrace{-\frac{1}{6},0,1,0,0,-1}_{\in\mathbb{R}^{6}},\ldots,\underbrace{-\frac{1}{6},0,1,0,0,-1}_{\in\mathbb{R}^{6}}\Big{]}^{T}\in\mathbb{R}^{3N}.$

    * –

$k=3$,
$u_{2}=\Big{[}\underbrace{0,-\frac{3}{20},0,1}_{\in\mathbb{R}^{4}},\underbrace{0,-\frac{3}{20},0,1}_{\in\mathbb{R}^{4}},\ldots,\underbrace{0,-\frac{3}{20},0,1}_{\in\mathbb{R}^{4}}\Big{]}^{T}\in\mathbb{R}^{4N}.$

    * –

$k=4$,
$v_{2}=\Big{[}\underbrace{-\frac{3}{280},0,\frac{3}{14},0,-1,0,0,-\frac{3}{14},0,1}_{\in\mathbb{R}^{10}},\ldots,\underbrace{-\frac{3}{280},0,\frac{3}{14},0,-1,0,0,-\frac{3}{14},0,1}_{\in\mathbb{R}^{10}}\Big{]}^{T}\in\mathbb{R}^{5N}.$

    * –

case $k=5$,
$u_{2}=\Big{[}\underbrace{0,\frac{5}{336},0,-\frac{5}{18},0,1}_{\in\mathbb{R}^{6}},\underbrace{0,\frac{5}{336},0,-\frac{5}{18},0,1}_{\in\mathbb{R}^{6}},\ldots,\underbrace{0,\frac{5}{336},0,-\frac{5}{18},0,1}_{\in\mathbb{R}^{6}}\Big{]}^{T}\in\mathbb{R}^{6N}.$

In the even case, when $k$ is even, we observe a double pattern which is
repeated $N/2$ (since when $k$ is even, $N$ is even to ensure $(k+1)N$ is
even). Once we get $(u_{1},u_{2})$, the formula $\Pi x=\langle x,u_{1}\rangle
u_{1}+\langle x,u_{2}\rangle u_{2}$ enables to get $\Pi$.

##  Appendix D Appendix D: Lawson-Fourier method for Vlasov-Maxwell 1dx-2dv.

In this appendix, we extend the method presented in [4] to the Vlasov-Maxwell
model in 1dx-2dv. This approach is compared to the Lawson-DG method in the
numerical section 5.

Starting from the Vlasov-Maxwell 1dx-2dv model (30) satisfied by
$f(t,x,v_{1},v_{2}),E_{1}(t,x),E_{2}(t,x),B(t,x)$, with $t\geq 0,x\in[0,L]$
and $(v_{1},v_{2})\in\mathbb{R}^{2}$, we shall use a Fourier method in the
space direction $x$ and we consider a grid in the velocity direction
$v_{\ell,j_{\ell}}=v_{\ell,\min}+j_{\ell}\Delta v_{\ell},\ell=1,2$. Denoting
$\hat{f}_{k,j_{1},j_{2}}(t)$ the spatial Fourier coefficients of
$f(t,x,v_{j_{1}},v_{j_{2}})$ and $(E_{1,k},E_{2,k},B_{k})(t)$ the Fourier
coefficients of $(E_{1},E_{2},B)(t,x)$ then gives the following semi-
discretized scheme for $k=1,\dots,N_{x}$ ($N_{x}$ being the number of grid
points in $[0,L]$)

$\left\\{\begin{aligned}
&\partial_{t}\hat{f}_{k,j_{1},j_{2}}+v_{1,j_{1}}ik\hat{f}_{k,j_{1},j_{2}}+\widehat{({\cal
F}{\cal D}f)}_{k,j_{1},j_{2}}=0,\;\;\mbox{ with }{\cal
F}=(E_{1}+Bv_{2},E_{2}-Bv_{1}),\\\
&\partial_{t}\hat{B}_{k}=-ik\hat{E}_{2,k},\\\
&\partial_{t}\hat{E}_{1,k}=-\sum_{j_{1},j_{2}}v_{1,j_{1}}\hat{f}_{k,j_{1},j_{2}}\Delta
v_{1}\Delta v_{2},\\\
&\partial_{t}\hat{E}_{2,k}=-ik\hat{B}_{k}-\sum_{j_{1},j_{2}}v_{2,j_{2}}\hat{f}_{k,j_{1},j_{2}}\Delta
v_{1}\Delta v_{2},\end{aligned}\right.$ (91)

with the initial condition
$\hat{f}_{k,j_{1},j_{2}},\hat{B}_{k}(0),\hat{E}_{1,k}(0),\hat{E}_{2,k}(0)$
satisfying the Poisson equation
$ik\hat{E}_{1,k}(0)=\sum_{j_{1},j_{2}}\hat{f}_{k,j_{1},j_{2}}(0)\Delta
v_{1}\Delta v_{2}$ for $k\neq 0$. Let denote $U=(\hat{\bf
f},\hat{B},\hat{E}_{2},\hat{E}_{1})\in\mathbb{M}_{N_{v_{1}}N_{v_{1}}+3,N_{v_{1}}N_{v_{1}}+3}(\mathbb{C})$,
then the previous system can be rewritten as

$\partial_{t}U=LU+N(U),$ (92)

with

$\displaystyle U$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{llll}\hat{f}_{k,\star,1}\\\
\hat{f}_{k,\star,2}\\\ \vdots\\\
\\!\\!\\!\hat{f}_{k,\star,N_{v_{2}}}\\!\\!\\!\\!\\\ \hat{B}_{k}\\\
\hat{E}_{2,k}\\\
\hat{E}_{1,k}\end{array}\right),\;\;\;\;\;\;N(U)=\left(\begin{array}[]{llllccc}-\widehat{({\cal
F}{\cal D}f)}_{k,\star,1}\\\ -\widehat{({\cal F}{\cal D}f)}_{k,\star,2}\\\
\vdots\\\ -\widehat{({\cal F}{\cal D}f)}_{k,\star,N_{v_{2}}}\\\ 0\\\ 0\\\
0\end{array}\right),\;\;\;\;\hat{f}_{k,\star,j_{2}}\in\mathbb{C}^{N_{v_{1}}},\begin{array}[]{llcc}\forall
k&=1,\dots,N_{x},\\\ \forall j_{2}&=1,\dots,N_{v_{2}},\end{array}$ (109)
$\displaystyle L\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\left(\begin{array}[]{lllllllccc}-ik\mbox{diag}(\vec{v_{1}})&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&{\bf 0}_{N_{v_{1}},N_{v_{1}}}\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&{\bf 0}_{N_{v_{1}},1}&\dots&{\bf 0}_{N_{v_{1}},1}\\\
{\bf 0}_{N_{v_{1}},N_{v_{1}}}&-ik\mbox{diag}(\vec{v_{1}})&{\bf
0}_{N_{v_{1}},N_{v_{1}}}\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\vdots&\vdots&\vdots\\\ \vdots&\ddots&\ddots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\vdots&\vdots&\vdots\\\ {\bf
0}_{N_{v_{1}},N_{v_{1}}}&\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&-ik\mbox{diag}(\vec{v_{1}})&{\bf
0}_{N_{v_{1}},1}&\vdots&{\bf 0}_{N_{v_{1}},1}\\\ {\bf
0}_{1,N_{v_{1}}}&\ldots&{\bf 0}_{1,N_{v_{1}}}&{\bf 0}_{1,N_{v_{1}}}&0&-ik&0\\\
-\Delta v_{1}\Delta v_{2}(\vec{v_{2}})_{1}{\bf 1}&-\Delta v_{1}\Delta
v_{2}(\vec{v_{2}})_{2}{\bf 1}&\ldots&-\Delta v_{1}\Delta
v_{2}(\vec{v_{2}})_{N_{vy}}{\bf 1}&-ik&0&0\\\ -\Delta v_{1}\Delta
v_{2}\vec{v_{1}}&-\Delta v_{1}\Delta v_{2}\vec{v_{1}}&\ldots&-\Delta
v_{1}\Delta v_{2}\vec{v_{1}}&0&0&0\end{array}\right)$ (117)

where we denote $\vec{v_{1}}$ the vector with components
$(\vec{v_{1}})_{j_{1}}=v_{1,\min}+j_{1}\Delta v_{1}$ and $\vec{v_{2}}$ the
vector with components $(\vec{v_{2}})_{j_{2}}=v_{2,\min}+j_{2}\Delta v_{2}$.
Moreover, diag$(\vec{v_{1}})$ denotes the diagonal matrix whose diagonal is
composed of $\vec{v_{1}}$, ${\bf 1}\in{\mathbb{R}}^{N_{v_{1}}}$ denotes the
vector with components $1$ and ${\bf 0}_{A,B}$ is a matrix with $A$ lines and
$B$ columns with zeros. The size of the matrix $L$ is $(N_{v_{1}}N_{v_{2}}+3)$
and in spite of its size, one can see that $L$ is sparse.

A key point is to compute $\exp(L\Delta t)$ to design an exponential scheme
approximating (92). Denoting $U^{n}\approx U(t^{n}),t^{n}=n\Delta t,\Delta
t>0$, the first order Lawson scheme is

$U^{n+1}=\exp(L\Delta t)U^{n}+\Delta t\exp(L\Delta t)N(U^{n}).$

The following proposition gives an explicit expression of $\exp(-\Delta tL)$.

######  Proposition 9.

The exponential of the matrix $L$ given by (117) is given by

$\left(\begin{array}[]{lllllllccc}e^{-ik\Delta t\mbox{diag}(\vec{v_{1}})}&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&{\bf 0}_{N_{v_{1}},N_{v_{1}}}\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&{\bf 0}_{N_{v_{1}},1}&{\bf 0}_{N_{v_{1}},1}&{\bf
0}_{N_{v_{1}},1}\\\ {\bf 0}_{N_{v_{1}},N_{v_{1}}}&e^{-ik\Delta
t\mbox{diag}(\vec{v_{1}})}&{\bf 0}_{N_{v_{1}},N_{v_{1}}}\ldots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\vdots&\vdots&\vdots\\\ \vdots&\ddots&\ddots&{\bf
0}_{N_{v_{1}},N_{v_{1}}}&\vdots&\vdots&\vdots\\\ {\bf
0}_{N_{v_{1}},N_{v_{1}}}&\ldots&{\bf 0}_{N_{v_{1}},N_{v_{1}}}&e^{-ik\Delta
t\mbox{diag}(\vec{v_{1}})}&{\bf 0}_{N_{v_{1}},1}&{\bf 0}_{N_{v_{1}},1}&{\bf
0}_{N_{v_{1}},1}\\\ {}_{e}{\cal B}_{k,1}&{}_{e}{\cal
B}_{k,2}&\ldots&{}_{e}{\cal B}_{k,N_{v_{2}}}&\cos(k\Delta t)&-i\sin(k\Delta
t)&0\\\ {}_{e}{\cal E}_{2,k,1}&{}_{e}{\cal E}_{2,k,2}&\ldots&{}_{e}{\cal
E}_{2,k,N_{v_{2}}}&-i\sin(k\Delta t)&\cos(k\Delta t)&0\\\ {}_{e}{\cal
E}_{1,k}&{}_{e}{\cal E}_{1,k}&\ldots&{}_{e}{\cal
E}_{1,k}&0&0&1\end{array}\right)$

where ${}_{e}{\cal B}_{k,j_{2}},{}_{e}{\cal E}_{2,k,j_{2}},{}_{e}{\cal
E}_{1,k}\in\mathbb{M}_{1,N_{v_{1}}}(\mathbb{C})$ for $j_{2}=1,\dots,N_{v_{2}}$
are given by

$\displaystyle{}_{e}{\cal B}_{k,j_{2}}$ $\displaystyle=$
$\displaystyle\frac{i\Delta v_{1}\Delta
v_{2}}{k}v_{2,j_{2}}\vec{\alpha},\;\;{}_{e}{\cal E}_{2,k,j_{2}}=\frac{i\Delta
v_{1}\Delta v_{2}}{k}v_{2,j_{2}}\vec{\beta},\;\;{}_{e}{\cal
E}_{1,k}=\frac{i\Delta v_{1}\Delta v_{2}}{k}({\bf 1}-e^{-ik\vec{v_{1}}}),$

where the vectors
$\vec{\alpha},\vec{\beta}\in\mathbb{M}_{1,N_{v_{1}}}(\mathbb{C})$ whose
components are given by

$\displaystyle{\alpha}_{j_{1}}$ $\displaystyle=$
$\displaystyle-\frac{e^{-ik\Delta t}}{2(1-v_{1,j_{1}})}-\frac{e^{ik\Delta
t}}{2(1+v_{1,j_{1}})}+\frac{e^{-ik\Delta tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})},$
$\displaystyle{\beta}_{j_{1}}$ $\displaystyle=$
$\displaystyle-\frac{e^{-ik\Delta t}}{2(1-v_{1,j_{1}})}+\frac{e^{ik\Delta
t}}{2(1+v_{1,j_{1}})}+\frac{v_{1,j_{1}}e^{-ik\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})}.$

###### Proof.

To compute $\exp(\Delta tL)$, we will solve exactly the linear part of the
Vlasov-Maxwell system. First of all, we observe that the components associated
to $\hat{f}$ is diagonal and can be solved independently so that the
$N_{v_{1}}\times N_{v_{2}}$ top left block of $\exp(\Delta tL)$ is diagonal
and is equal to $e^{-ik\Delta t\vec{v_{1}}}$ ($N_{v_{2}}$ times). Second, the
$3\times 3$ right bottom block corresponds to the homogeneous Maxwell
equations and its exponential can be computed as

$\exp\left(\begin{array}[]{lll}0&-ik\Delta t&0\\\ -ik\Delta t&0&0\\\
0&0&0\end{array}\right)=\left(\begin{array}[]{lll}\cos(k\Delta
t)&i\sin(k\Delta t)&0\\\ i\sin(k\Delta t)&\cos(k\Delta t)&0\\\
0&0&1\end{array}\right).$

It remains to compute the three last lines corresponding to the coupling
between the Vlasov and Maxwell equations.

##### Computation of ${}_{e}{\cal E}_{1,k}$: solve $\hat{E}_{1,k}$

First, we have for $\hat{f}_{k,j_{1},j_{2}}(t)$

$\hat{f}_{k,j_{1},j_{2}}(t)=e^{-ikv_{1,j_{1}}(t-t^{n})}\hat{f}_{k,j_{1},j_{2}}(t^{n}),$

which enables to compute $\hat{E}_{1,k}(t^{n+1})$

$\displaystyle\hat{E}_{1,k}(t^{n+1})$ $\displaystyle=$
$\displaystyle\hat{E}_{1,k}(t^{n})-\sum_{j_{1},j_{2}}\int_{t^{n}}^{t^{n+1}}e^{-ikv_{1,j_{1}}(t-t^{n})}dtv_{1,j_{1}}{f}_{k,j_{1},j_{2}}(t^{n})\Delta
v_{1}\Delta v_{2}$ $\displaystyle=$
$\displaystyle\hat{E}_{1,k}(t^{n})+\frac{i\Delta v_{1}\Delta
v_{2}}{k}\sum_{j_{1},j_{2}}(1-e^{-ikv_{1,j_{1}}\Delta
t})\hat{f}_{k,j_{1},j_{2}}(t^{n}),$

from which we deduce the expression of ${}_{e}{\cal E}_{1,k}$.

##### Computation of ${}_{e}{\cal B}_{k},{}_{e}{\cal E}_{2,k}$: solve
$\hat{B}_{k},\hat{E}_{2,k}$

Next, we focus on the calculation of $\hat{E}_{2,k}(t^{n+1})$ and
$\hat{B}_{k}(t^{n+1})$. We write down the equations

$\displaystyle\frac{d}{dt}\hat{E}_{2,k}(t)$ $\displaystyle=$ $\displaystyle-
ik\hat{B}_{k}(t)-\sum_{j_{1},j_{2}}e^{-ikv_{1,j_{1}}(t-t^{n})}v_{2,j_{1}}{f}_{k,j_{1},j_{2}}(t^{n})\Delta
v_{1}\Delta v_{2}$ $\displaystyle\frac{d}{dt}\hat{B}_{k}(t)$ $\displaystyle=$
$\displaystyle-ik\hat{E}_{2,k}(t)$

which can be rewritten as $\frac{dU}{dt}=MU+R$ with
$U(t)=(\hat{E}_{2,k}(t),\hat{B}_{k}(t))$ and

$M=\left(\begin{array}[]{ll}0&-ik\\\ -ik&0\end{array}\right)\mbox{ and
}R=\left(\begin{array}[]{ll}R_{1}(t)\\\ 0\end{array}\right),\mbox{ with
}R_{1}(t)=-\sum_{j_{1},j_{2}}e^{-ikv_{1,j_{1}}(t-t^{n})}v_{2,j_{1}}\hat{f}_{k,j_{1},j_{2}}(t^{n})\Delta
v_{1}\Delta v_{2}.$

Thus, one can write the variation of constant formula

$U(t^{n+1})=e^{M\Delta
t}U(t^{n})+\int_{t^{n}}^{t^{n+1}}e^{-M(t-t^{n+1})}R(t)dt.$ (118)

First, $e^{M\Delta t}$ reads as

$e^{M\Delta t}=\exp\left(\begin{array}[]{ll}0&-ik\Delta t\\\ -ik\Delta
t&0\end{array}\right)=\left(\begin{array}[]{ll}\cos(kt)&-i\sin(kt)\\\
-i\sin(kt)&\cos(kt)\end{array}\right).$

Second, one has to compute the integral term in (118)

$\displaystyle\int_{t^{n}}^{t^{n+1}}e^{-M(t-t^{n+1})}R(t)dt$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ll}\int_{t^{n}}^{t^{n+1}}\cos(k(t-t^{n+1}))R_{1}(t)dt\\\
\int_{t^{n}}^{t^{n+1}}i\sin(k(t-t^{n+1}))R_{1}(t)dt\\\ \end{array}\right)$
(121) $\displaystyle=$
$\displaystyle-\left(\begin{array}[]{ll}\sum_{j_{1},j_{2}}{\cal
I}_{1}v_{2,j_{1}}\hat{f}^{n}_{k,j_{1},j_{2}}\Delta v_{1}\Delta v_{2}\\\
\sum_{j_{1},j_{2}}{\cal I}_{2}v_{2,j_{1}}\hat{f}^{n}_{k,j_{1},j_{2}}\Delta
v_{1}\Delta v_{2}\end{array}\right)$ (124)

where ${\cal I}_{1}$ and ${\cal I}_{2}$ are given by

$\displaystyle{\cal I}_{1}$ $\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\Big{[}\cos(k(t-t^{n+1}))e^{-ikv_{1,j_{1}}(t-t^{n})}\Big{]}dt=\frac{ie^{-ik\Delta
t}}{2k(1-v_{1,j_{1}})}-\frac{ie^{ik\Delta
t}}{2k(1+v_{1,j_{1}})}-\frac{iv_{1,j_{1}}e^{-ik\Delta
tv_{1,j_{1}}}}{k(1-v_{1,j_{1}}^{2})},$ $\displaystyle{\cal I}_{2}$
$\displaystyle=$
$\displaystyle\int_{t^{n}}^{t^{n+1}}\Big{[}i\sin(k(t-t^{n+1}))e^{-ikv_{1,j_{1}}(t-t^{n})}\Big{]}dt=\frac{ie^{-ik\Delta
t}}{2k(1-v_{1,j_{1}})}+\frac{ie^{ik\Delta
t}}{2k(1+v_{1,j_{1}})}-\frac{ie^{-ik\Delta
tv_{1,j_{1}}}}{k(1-v_{1,j_{1}}^{2})}.$

Inserting these calculations in (118) leads to the following expression for
$\hat{E}_{2,k}(t)$

$\displaystyle\hat{E}_{2,k}(t^{n+1})$ $\displaystyle=$
$\displaystyle\cos(k\Delta t)\hat{E}_{2,k}(t^{n})-i\sin(k\Delta
t)\hat{B}_{k}(t^{n})+\frac{i\Delta v_{1}\Delta
v_{2}}{k}\sum_{j_{1},j_{2}}v_{2,j_{2}}\beta_{j_{1}}\hat{f}_{k,j_{1},j_{2}}(t^{n}),$
$\displaystyle\hat{B}_{k}(t^{n+1})$ $\displaystyle=$
$\displaystyle-i\sin(k\Delta t)E_{2,k}(t^{n})+\cos(k\Delta
t)B_{k}(t^{n})+\frac{i\Delta v_{1}\Delta
v_{2}}{k}\sum_{j_{1},j_{2}}v_{2,j_{2}}\alpha_{j_{1}}\hat{f}_{k,j_{1},j_{2}}(t^{n}),$

where
$\vec{\beta}=[\beta_{1},\beta_{2},\dots,\beta_{N_{v_{1}}}]\in\mathbb{C}^{N_{v_{1}}}$
and
$\vec{\alpha}=[\alpha_{1},\alpha_{2},\dots,\alpha_{N_{v_{1}}}]\in\mathbb{C}^{N_{v_{1}}}$
are given by

$\displaystyle\beta_{j_{1}}$ $\displaystyle=$
$\displaystyle-\frac{e^{-ik\Delta t}}{2(1-v_{1,j_{1}})}+\frac{e^{ik\Delta
t}}{2(1+v_{1,j_{1}})}+\frac{v_{1,j_{1}}e^{-ik\Delta
tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})},$ $\displaystyle\alpha_{j_{1}}$
$\displaystyle=$ $\displaystyle-\frac{e^{-ik\Delta
t}}{2(1-v_{1,j_{1}})}-\frac{e^{ik\Delta
t}}{2(1+v_{1,j_{1}})}+\frac{e^{-ik\Delta tv_{1,j_{1}}}}{(1-v_{1,j_{1}}^{2})}.$

We conclude by writing the vectors ${}_{e}{\cal B}_{k},{}_{e}{\cal E}_{2,k}$
corresponding to $\hat{B}_{k}$ and $\hat{E}_{2,k}$

${}_{e}{\cal B}_{k,j_{2}}=\frac{i\Delta v_{1}\Delta
v_{2}}{k}v_{2,j_{2}}\vec{\alpha},\;\;{}_{e}{\cal E}_{2,k,j_{2}}=\frac{i\Delta
v_{1}\Delta v_{2}}{k}v_{2,j_{2}}\vec{\beta}.$

∎

## References

  * [1] J. W. Banks and J. A. F. Hittinger.  A new class of nonlinear finite-volume methods for vlasov simulation.  IEEE Transactions on Plasma Science, 38(9):2198–2207, 2010.
  * [2] J. W. Banks, A. G. Odu, R. Berger, T. Chapman, W. Arrighi, and S. Brunner.  High-order accurate conservative finite difference methods for vlasov equations in 2d+ 2v.  SIAM Journal on Scientific Computing, 41(5):B953–B982, 2019.
  * [3] R. Biswas, K. D. Devine, and J. E. Flaherty.  Parallel, adaptive finite element methods for conservation laws.  Applied Numerical Mathematics, 14(1-3):255–283, 1994.
  * [4] B. Boutin, A. Crestetto, N. Crouseilles, and J. Massot.  Modified lawson methods for vlasov equations.  2022\.
  * [5] X. Cai, J.-M. Qiu, and Y. Yang.  An eulerian-lagrangian discontinuous galerkin method for transport problems and its application to nonlinear dynamics.  Journal of Computational Physics, 439:110392, 2021.
  * [6] M. Caliari, F. Cassini, L. Einkemmer, A. Ostermann, and F. Zivcovich.  A $\mu$-mode integrator for solving evolution equations in kronecker form.  Journal of Computational Physics, 455:110989, 2022.
  * [7] F. Califano, F. Pegoraro, and S. V. Bulanov.  Spatial structure and time evolution of the weibel instability in collisionless inhomogeneous plasmas.  Physical review E, 56(1):963, 1997.
  * [8] M. A. Celia, T. F. Russell, I. Herrera, and R. E. Ewing.  An eulerian-lagrangian localized adjoint method for the advection-diffusion equation.  Advances in water resources, 13(4):187–206, 1990.
  * [9] Y. Cheng, I. M. Gamba, F. Li, and P. J. Morrison.  Discontinuous galerkin methods for the vlasov–maxwell equations.  SIAM Journal on Numerical Analysis, 52(2):1017–1049, 2014.
  * [10] P. G. Ciarlet.  The finite element method for elliptic problems.  SIAM, 2002.
  * [11] B. Cockburn and C.-W. Shu.  Runge–kutta discontinuous galerkin methods for convection-dominated problems.  Journal of scientific computing, 16:173–261, 2001.
  * [12] A. Crestetto, N. Crouseilles, Y. Li, and J. Massot.  Comparison of high-order eulerian methods for electron hybrid model.  Journal of Computational Physics, 451:110857, 2022.
  * [13] M. Croci and J. Muñoz-Matute.  Exploiting kronecker structure in exponential integrators: Fast approximation of the action of $\varphi$-functions of matrices via quadrature.  Journal of Computational Science, 67:101966, 2023.
  * [14] N. Crouseilles, L. Einkemmer, and J. Massot.  Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations.  Journal of Computational Physics, 420:109688, 2020.
  * [15] L. Einkemmer.  A performance comparison of semi-lagrangian discontinuous galerkin and spline based vlasov solvers in four dimensions.  Journal of Computational Physics, 376:937–951, 2019.
  * [16] F. Filbet and E. Sonnendrücker.  Comparison of eulerian vlasov solvers.  Computer Physics Communications, 150(3):247–266, 2003.
  * [17] M. Hochbruck and A. Ostermann.  Exponential integrators.  Acta Numerica, 19:209–286, 2010.
  * [18] X. Hong and J.-M. Qiu.  A generalized eulerian-lagrangian discontinuous galerkin method for transport problems.  Journal of Computational Physics, 464:111160, 2022.
  * [19] C. Johnson and J. Pitkäranta.  An analysis of the discontinuous galerkin method for a scalar hyperbolic equation.  Mathematics of computation, 46(173):1–26, 1986.
  * [20] A. J. Klimas and W. M. Farrell.  A splitting algorithm for vlasov simulation with filamentation filtration.  Journal of computational physics, 110(1):150–163, 1994.
  * [21] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker.  Gempic: geometric electromagnetic particle-in-cell methods.  Journal of Plasma Physics, 83(4):905830401, 2017.
  * [22] P. Lesaint and P.-A. Raviart.  On a finite element method for solving the neutron transport equation.  Publications mathématiques et informatique de Rennes, (S4):1–40, 1974.
  * [23] Y. Liu, C.-W. Shu, and M. Zhang.  Sub-optimal convergence of discontinuous galerkin methods with central fluxes for linear hyperbolic equations with even degree polynomial approximations.  Journal of Computational Mathematics, 39:629–648, 2021.
  * [24] J. Muñoz-Matute, D. Pardo, and V. M. Calo.  Exploiting the kronecker product structure of $\varphi$\- functions in exponential integrators.  International Journal for Numerical Methods in Engineering, 123(9):2142–2161, 2022.
  * [25] T. E. Peterson.  A note on the convergence of the discontinuous galerkin method for a scalar hyperbolic equation.  SIAM Journal on Numerical Analysis, 28(1):133–140, 1991.