# Exponential DG methods for Vlasov equations

Nicolas Crouseilles, Xue Hong

► To cite this version:

Nicolas Crouseilles, Xue Hong. Exponential DG methods for Vlasov equations. Journal of Computational Physics, 2024, 498, pp.1-38. 10.1016/j.jcp.2023.112682. hal-04389308

HAL Id: hal-04389308

https://hal.science/hal-04389308v1

Submitted on 11 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Exponential DG methods for Vlasov equations

Nicolas Crouseilles¹ and Xue Hong²

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

²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.

**Key words:** exponential Lawson Runge-Kutta method; discontinuous Galerkin method; Vlasov-Poisson equation; Vlasov Maxwell equations; conservative property.

## 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 of 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 \ge 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 [35, 26] 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 [25, 21] or finite differences/volumes methods [3, 2]. 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 [31, 33, 20] or
arbitrary Lagrangian-Eulerian methods [9, 36, 6, 23], 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 [16]).

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 ([22] 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 [18, 16, 5], 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 [15]) which is important for physically relevant problems (applications of the DG method to kinetic models also exist, e.g. [14, 12, 1, 19, 10, 11]); moreover, one element only communicates with its immediate neighbors which is very important for parallelization capability [4] 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:

$$ \begin{cases} u_t + au_x = 0, \quad x \in [x_a, x_b], \\ u(0, x) = u_0(x). \end{cases} \tag{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, \quad j = 1, \dots, N, \tag{2} $$

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

$$ \text{central flux: } \hat{F}|_{x_{j \pm \frac{1}{2}}} = \frac{u^{-} + u^{+}}{2} |_{x_{j \pm \frac{1}{2}}}, \quad \text{upwind flux: } \hat{F}|_{x_{j \pm \frac{1}{2}}} = u^{-}|_{x_{j \pm \frac{1}{2}}}, \tag{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_j^m(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_j^m(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), \tag{4} $$

the semi-discretized problem simply becomes

$$ \frac{d\mathbf{u}}{dt} = A\mathbf{u}, \quad \text{with } A \in \mathbb{M}_{(k+1)N, (k+1)N}(\mathbb{R}). \tag{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 [15, 38, 39]. 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, \quad \forall n \in \mathbb{N}, \tag{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_h^n \in V_h^k$ from

$$u_h^n(x) = \sum_{m=0}^{k} (u_j^m)^n \xi_j^m(x), \quad \forall x \in I_j, \quad j = 1, \dots, N. \qquad (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} \begin{pmatrix} C_1 & C_2 & \mathbf{0} & \dots & C_3 \\ C_3 & C_1 & C_2 & \mathbf{0} & \dots \\ \mathbf{0} & \ddots & \ddots & \ddots & C_2 \\ C_2 & \mathbf{0} & \dots & C_3 & C_1 \end{pmatrix}, \quad \text{with } C_j = M^{-1}D_j \quad (j = 1, 2, 3) \qquad (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)

$$M_{\ell,m} = \frac{(1/2)^{m+\ell-1}}{m+\ell-1} [1 - (-1)^{m+\ell-1}], \quad (D_2)_{\ell,m} = (-1)^m (1/2)^{m+\ell-1}, \\ (D_1)_{\ell,m} = (1/2)^{m+\ell-2} \left( \frac{\ell-1}{m+\ell-2} - \frac{1}{2} \right) [1 - (-1)^{m+\ell-2}], \quad (D_3)_{\ell,m} = (-1)^{\ell-1} (1/2)^{m+\ell-1}, \qquad (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 [34], 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 $C_j = C_1 + \rho_j C_2 + \rho_j^{N-1} C_3$ for $j = 0, \dots, N-1$, where $C_j \in \mathbb{M}_{(k+1),(k+1)}$. Then, we checked numerically that the eigenvalues of $C_j$ are pure imaginary for all $j = 0, \dots, N-1$, and we deduce from [34] 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 $\{0, \lambda_j, \bar{\lambda}_j, j = 1, \dots, d\}$, in particular $0$ is a simple eigenvalue in this case. Since $P_A(-\lambda) = P_A(\lambda)$, we deduce $\mathrm{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 $\{0, \lambda_j, \bar{\lambda}_j, j = 1, \dots, d-1\}$, in particular $0$ is a double eigenvalue in this case. Since $P_A(\lambda) = P_A(-\lambda)$, we deduce $\mathrm{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 us 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)\| \le C$, with $\|\cdot\|$ an induced matrix norm.

*Proof.* First, we write $A = \Delta x^{-1}A_1$ from (8) with $A_1$ independent of $\Delta x$ (with $\Delta x = (x_b - x_a)/N$) and since $A$ is diagonalizable from Remarks 1 & 2, 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 does not depend on $N$), so does the matrix $P$. Thus, there exists $C > 0$ (independent of $N$) such that $\text{cond}(P) = \|P\|\|P^{-1}\| \le 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}\| \le \text{cond}(P) \|\exp(\Delta x^{-1}D)\| = \text{cond}(P) \le 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} \le 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} \le C\Delta x^k$, where $C$ is a positive constant independent on $\Delta x$ (see the details in [32, 27, 24, 30, 38, 39], 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} \le \|u(t^n, \cdot) - u_h(t^n, \cdot)\|_{L^2} + \|u_h(t^n, \cdot) - u_h^n\|_{L^2} \le 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 [24, 30, 38, 39]). However, $k$th order sub-optimal convergence rate is proved for DG with central flux in [28] 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) \ge 0$ and the electric field $E(t, x) \in \mathbb{R}$ with $t \ge 0, x \in [0, L] (L > 0)$ and $v \in \mathbb{R}$,

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial f}{\partial t} + v \frac{\partial f}{\partial x} + E \frac{\partial f}{\partial v} = 0, \\[1.5em] \displaystyle \partial_t E = - \int_{\mathbb{R}} v f dv + \bar{J}, \quad \bar{J} = \frac{1}{L} \int_0^L \int_{\mathbb{R}} v f dx dv, \end{array} \right. \qquad (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 dv dx$ and periodic boundary conditions are imposed in space. The Vlasov-Ampère system is equivalent to the Vlasov-Poisson model where the electric field satisfies the Poisson equation $\partial_x E = \int_{\mathbb{R}} f dv - \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-1/2}, x_{i+1/2}]$, $i = 1, \dots, 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), \quad j = 0, \dots, N_v. $$

As in the previous section, we denote $f_j(t) \in \mathbb{R}^{(k+1)N_x}$ the vector of the DG coefficients $f_i^m(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_i^m$ 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$:

$$ \sum_{m=0}^{k} \left[ \left( \partial_t f_i^m(t, v_j) \xi^m, \xi^\ell \right)_{I_i} - v_j \left( f_i^m(t, v_j) \xi^m, \partial_x \xi^\ell(x) \right)_{I_i} \right] + v_j \left[ \{f_h(t, x, v_j)\} \xi^\ell \right]_{i-1/2}^{i+1/2} \qquad (11) \\ + \sum_{m=0}^{k} \left( \sum_{n=0}^{k} E_i^n \xi^n(\mathcal{D} f_i^m)(t, v_j) \xi^m, \xi^\ell \right)_{I_i} = 0, $$

where we used the central flux $\{f_h(t, x, v_j)\}|_{x_{i+1/2}} = \frac{1}{2}(f_h(t, x_{i+1/2}^+, v_j) + f_h(t, x_{i+1/2}^-, v_j))$, $\ell = 0, 1, 2, \dots, k$, $i = 1, 2, \dots, 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), \dots, f_i^k(t, v_j))^T, \quad i = 1, 2, \dots, N_x, $$
$$ \mathbf{E}_i(t) = (E_i^0(t), E_i^1(t), \dots, E_i^k(t))^T, \quad i = 1, 2, \dots, N_x, $$

and

$$ \mathcal{D}\mathbf{f}_{j,i}(t) = ((\mathcal{D}f_i^0)(t, v_j), (\mathcal{D}f_i^1)(t, v_j), \dots, (\mathcal{D}f_i^k)(t, v_j))^T, \quad i = 1, 2, \dots, 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$.

$$ \Delta x \begin{pmatrix} M & & & \\ & M & & \\ & & \ddots & \\ & & & M \end{pmatrix} \frac{d}{dt} \begin{pmatrix} \mathbf{f}_{j,1} \\ \mathbf{f}_{j,2} \\ \vdots \\ \mathbf{f}_{j,N_x} \end{pmatrix} - v_j \begin{pmatrix} D_1 & D_2 & \dots & D_3 \\ D_3 & D_1 & D_2 & \\ & \ddots & \ddots & D_2 \\ D_2 & \dots & D_3 & D_1 \end{pmatrix} \begin{pmatrix} \mathbf{f}_{j,1} \\ \mathbf{f}_{j,2} \\ \vdots \\ \mathbf{f}_{j,N_x} \end{pmatrix} + \begin{pmatrix} B_1 & & & \\ & B_2 & & \\ & & \ddots & \\ & & & B_{N_x} \end{pmatrix} \begin{pmatrix} (\mathcal{D}\mathbf{f})_{j,1} \\ (\mathcal{D}\mathbf{f})_{j,2} \\ \vdots \\ (\mathcal{D}\mathbf{f})_{j,N_x} \end{pmatrix} = \mathbf{0}, \quad (12) $$

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$

$$ \mathbf{f}_j(t) = (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, $$

$$ \mathbf{E}(t) = (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 $\mathcal{D}\mathbf{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}}(\mathcal{D}\mathbf{f})_j, \quad (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 x M)^{-1} B_i, i = 1, \dots, N_x$. Let us remark that we consider central finite differences method to approximate $(\mathcal{D}\mathbf{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)dv - \bar{\rho}$ would be $-A\mathbf{E}^0 = \sum_j \mathbf{f}_j^0 \Delta v - \bar{\rho}$ (with $\mathbf{E}^0$ and $\mathbf{f}_j^0$ 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)dx = 0$ for Poisson equation itself. We then consider the following discretized Poisson equation

$$ -(A + \Pi)\mathbf{E}^0 = (\mathbf{1} - \Pi)\left(\sum_j \mathbf{f}_j^0 \Delta v - \bar{\rho}\right) = (\mathbf{1} - \Pi)\sum_j \mathbf{f}_j^0 \Delta v, \quad (14) $$

with $\mathbf{1}$ the identity matrix of size $(k+1)N_x$ and where in the last equality, we used the fact that constants belong to $\text{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) = (\mathbf{1} - \Pi) \sum_j \mathbf{f}_j(t) \Delta v$. Considering
the time derivative of the latter equation and using (13) leads to

$$ -(A + \Pi)\partial_t \mathbf{E}(t) = (\mathbf{1} - \Pi)\partial_t \sum_j \mathbf{f}_j(t)\Delta v = (\mathbf{1} - \Pi)\sum_j v_j A \mathbf{f}_j(t)\Delta v \\ = A(\mathbf{1} - \Pi)\sum_j v_j \mathbf{f}_j(t)\Delta v = (A + \Pi)(\mathbf{1} - \Pi)\sum_j v_j \mathbf{f}_j(t)\Delta v, $$

where we used $\sum_j \mathcal{D}\mathbf{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) = -(\mathbf{1} - \Pi) \sum_j v_j \mathbf{f}_j(t) \Delta v. \qquad (15) $$

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

$$ \left\{ \begin{array}{l} \partial_t \mathbf{f}_j = v_j A \mathbf{f}_j - \tilde{\mathbf{E}}(\mathcal{D} \mathbf{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{array} \right. \qquad (16) $$

In view of the time discretization, we introduce the following vector of semi-discrete unknown $U = (\tilde{\mathbf{f}}, \mathbf{E}) \in \mathbb{R}^{(k+1)N_x(N_v+1)}$ with $\tilde{\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), \qquad (17) $$

with $L \in \mathbb{M}_{(k+1)N_x(N_v+1),(k+1)N_x(N_v+1)}$ given by

$$ L = \begin{pmatrix} v_1 A & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} & \dots & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} \\ \mathbf{0}_{\tilde{N},\tilde{N}} & v_2 A & \mathbf{0}_{\tilde{N},\tilde{N}} & \dots & \mathbf{0}_{\tilde{N},\tilde{N}} & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ \mathbf{0}_{\tilde{N},\tilde{N}} & \dots & \mathbf{0}_{\tilde{N},\tilde{N}} & v_{N_v} A & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} \\ -\Delta v v_1 (\mathbf{1} - \Pi) & -\Delta v v_2 (\mathbf{1} - \Pi) & \dots & -\Delta v v_{N_v} (\mathbf{1} - \Pi) & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} \end{pmatrix} \qquad (18) $$

where we denote $\tilde{N} = (k+1)N_x$ and $\mathbf{1} = \mathbf{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

$$ U = (\mathbf{f}_1, \mathbf{f}_2, \dots, \mathbf{f}_{N_v}, \mathbf{E})^T, \qquad N(U) = (-\tilde{\mathbf{E}}(\mathcal{D}\mathbf{f})_1, -\tilde{\mathbf{E}}(\mathcal{D}\mathbf{f})_2, \dots, -\tilde{\mathbf{E}}(\mathcal{D}\mathbf{f})_{N_v}, \mathbf{0})^T. $$

### 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 [22, 18] 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 [18].

Denoting $U^n = (\tilde{\mathbf{f}}^n, \mathbf{E}^n) \approx (\tilde{\mathbf{f}}(t^n), \mathbf{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 t L) U^n + \Delta t \exp(\Delta t L) N(U^n). \qquad (19) $$
High order methods can be obtained from Runge-Kutta methods using the corresponding Butcher tableau [18].

The key point of exponential methods lies in the computation of $\exp(\Delta tL)$ with $L$ given by (18). To compute $\exp(At) \forall t > 0$, one considers the linear part only $\partial_t U = LU$. First, we observe that the distribution function $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

$$f_j(t) = \exp(v_j At)f_j(0). \qquad (20)$$

Now, let us consider the electric field equation

$$\partial_t \mathbf{E} = -\sum_j v_j(\mathbf{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(\mathbf{1} - \Pi)\mathbf{f}_j(s)\Delta v ds.$$

Now, replacing $f_j(s)$ by $\exp(v_jAs)f_j(0)$ from (20) enables to get an explicit expression of $\mathbf{E}(t)$. Indeed, using properties used to derive (14) we have

$$
\begin{aligned}
\mathbf{E}(t) &= \mathbf{E}(0) - \sum_j \int_0^t v_j(\mathbf{1} - \Pi)\mathbf{f}_j(s)\Delta v ds = \mathbf{E}(0) - \sum_j \int_0^t v_j(\mathbf{1} - \Pi)\exp(v_jAs)\mathbf{f}_j(0)ds\Delta v \\
&= \mathbf{E}(0) - \sum_j \int_0^t v_j(\mathbf{1} - \Pi) \left[ \sum_{k=0}^{\infty} \frac{(v_jAs)^k}{k!} (\Pi\mathbf{f}_j(0) + (\mathbf{1} - \Pi)\mathbf{f}_j(0)) \right] ds\Delta v \\
&= \mathbf{E}(0) - \sum_j v_j(\mathbf{1} - \Pi) \int_0^t \left[ \Pi\mathbf{f}_j(0) + \sum_{k=0}^{\infty} \frac{(v_j(A + \Pi)s)^k}{k!} (\mathbf{1} - \Pi)\mathbf{f}_j(0) \right] ds\Delta v \\
&= \mathbf{E}(0) - \sum_j v_j(\mathbf{1} - \Pi) \left[ t\Pi\mathbf{f}_j(0) + \sum_{k=0}^{\infty} \frac{(v_j(A + \Pi))^k t^{k+1}}{(k+1)!} (\mathbf{1} - \Pi)\mathbf{f}_j(0) \right] \Delta v \\
&= \mathbf{E}(0) - \sum_j \left[ (\mathbf{1} - \Pi) \sum_{k=0}^{\infty} (A + \Pi)^{-1} \frac{(v_j(A + \Pi)t)^{k+1}}{(k+1)!} (\mathbf{1} - \Pi)\mathbf{f}_j(0) \right] \Delta v \\
&= \mathbf{E}(0) - \sum_j (\mathbf{1} - \Pi)(A + \Pi)^{-1}(\exp(v_jAt) - \mathbf{1})(\mathbf{1} - \Pi)\mathbf{f}_j(0)\Delta v \\
&= \mathbf{E}(0) - \sum_j (A + \Pi)^{-1}(\exp(v_jAt) - \mathbf{1})(\mathbf{1} - \Pi)\mathbf{f}_j(0)\Delta v \\
&= \mathbf{E}(0) + \Delta v \sum_j (A + \Pi)^{-1}(\mathbf{1} - \exp(v_jAt))\mathbf{f}_j(0).
\end{aligned}
$$

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

$$\left( \begin{array}{ccccccc}
e^{v_1 At} & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} & \dots & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} \\
\mathbf{0}_{\tilde{N},\tilde{N}} & e^{v_2 At} & \mathbf{0}_{\tilde{N},\tilde{N}} & \dots & \mathbf{0}_{\tilde{N},\tilde{N}} & \vdots & \vdots \\
\vdots & \ddots & \ddots & \dots & \vdots & \vdots & \vdots \\
\mathbf{0}_{\tilde{N},\tilde{N}} & \dots & \mathbf{0}_{\tilde{N},\tilde{N}} & \dots & e^{v_{Nv} At} & \mathbf{0}_{\tilde{N},\tilde{N}} & \mathbf{0}_{\tilde{N},\tilde{N}} \\
\Delta v \tilde{A}^{-1}(\mathbf{1} - e^{v_1 At}) & \Delta v \tilde{A}^{-1}(\mathbf{1} - e^{v_2 A\Delta t}) & \dots & \dots & \Delta v \tilde{A}^{-1}(\mathbf{1} - e^{v_{Nv} At}) & \mathbf{1} & \mathbf{1}
\end{array} \right). \qquad (21)$$
Hence, the exponential-DG scheme for the Vlasov-Ampère equation corresponds to (19)-(21). 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 (19)-(21) (and its generalization to high order Lawson Runge-Kutta) satisfied by $U^n = (\mathbf{f}^n, \mathbf{E}^n)$ preserves the following discretized Poisson equation:

$$ (A + \Pi)\mathbf{E}^n = -\sum_j (\mathbf{1} - \Pi)\mathbf{f}_j^n \Delta v, \quad \forall n \in \mathbb{N}^*, $$

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 $\text{Ker}(A)$ and $\mathbf{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 (\mathbf{1} - \Pi)\mathbf{f}_j^0 \Delta v$ holds at the initial time.

Next, from the scheme (19) with $\exp(Lt)$ given by (21), 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{\mathcal{E}} \mathcal{D} \mathbf{f})_j^n, \quad (22) $$

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

$$ \begin{aligned} \mathbf{E}^{n+1} &= \mathbf{E}^n + \Delta v \sum_j (A + \Pi)^{-1} (\mathbf{1} - \exp(v_j A \Delta t)) \mathbf{f}_j^n \\ &\quad - \Delta t \Delta v \sum_j (A + \Pi)^{-1} (\mathbf{1} - \exp(v_j A \Delta t)) (\tilde{\mathcal{E}} \mathcal{D} \mathbf{f})_j^n \\ &= \mathbf{E}^n + \Delta v \sum_j (A + \Pi)^{-1} (\mathbf{1} - \Pi) (\mathbf{1} - \exp(v_j A \Delta t)) \mathbf{f}_j^n \\ &\quad - \Delta t \Delta v \sum_j (A + \Pi)^{-1} (\mathbf{1} - \Pi) (\mathbf{1} - \exp(v_j A \Delta t)) (\tilde{\mathcal{E}} \mathcal{D} \mathbf{f})_j^n. \end{aligned} $$

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

$$ \begin{aligned} \Delta t \Delta v (A + \Pi)^{-1} \sum_j (\mathbf{1} - \Pi) \exp(v_j A \Delta t) (\tilde{\mathcal{E}} \mathcal{D} \mathbf{f})_j^n &= \Delta v (A + \Pi)^{-1} \sum_j (\mathbf{1} - \Pi) (\exp(v_j A \Delta t) \mathbf{f}_j^n - \mathbf{f}_j^{n+1}) \\ &= \Delta v (A + \Pi)^{-1} \sum_j \left[ (\mathbf{1} - \Pi) \exp(v_j A \Delta t) \mathbf{f}_j^n - (\mathbf{1} - \Pi) \mathbf{f}_j^{n+1} \right]. \end{aligned} $$

Finally, we then get

$$ \begin{aligned} \mathbf{E}^{n+1} &= \mathbf{E}^n + \Delta v \sum_j (A + \Pi)^{-1} (\mathbf{1} - \Pi) (\mathbf{1} - \exp(v_j A \Delta t)) \mathbf{f}_j^n \\ &\quad + \Delta v (A + \Pi)^{-1} \sum_j \left[ (\mathbf{1} - \Pi) \exp(v_j A \Delta t) \mathbf{f}_j^n - (\mathbf{1} - \Pi) \mathbf{f}_j^{n+1} \right] \\ &= \mathbf{E}^n + \Delta v \sum_j (A + \Pi)^{-1} ((\mathbf{1} - \Pi) \mathbf{f}_j^n - (\mathbf{1} - \Pi) \mathbf{f}_j^{n+1}). \end{aligned} \quad (23) $$

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.

<figure><img src="image_2.png" /></figure>
**Remark 4.** In this remark, we discuss how to compute the matrix Π. 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 Π 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 \ge 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 \, dv_1 \, dv_2 + \bar{J}_1, \quad \partial_t E_2 = -\partial_x B - \int_{\mathbb{R}^2} v_2 f \, dx \, dv_1 \, dv_2 + \bar{J}_2, \end{aligned} \right. \qquad (24) $$

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) dx dv_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 f \, dx \, dv_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}^*$. The definitions of the different objects are a direct extension of the definitions introduced in the previous part. Indeed, we denote by $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{\mathcal{F}}\mathcal{D}f)_{j_1,j_2}$ (with $\tilde{\mathcal{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 $(\mathcal{F} \cdot \nabla v)f$, with $\mathcal{F} = (E_1 + Bv_2, E_2 - Bv_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{\mathcal{F}} \mathcal{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} (\mathbf{1} - \Pi) \mathbf{f}_{j_1,j_2} \Delta v_1 \Delta v_2, \quad \partial_t \mathbf{E}_2 = A \mathbf{B} - \sum_{j_1,j_2} v_{2,j_2} (\mathbf{1} - \Pi) \mathbf{f}_{j_1,j_2} \Delta v_1 \Delta v_2, \end{aligned} \right. \qquad (25) $$

where A is given by (8)-(9) and **1** the identity matrix of size $(k+1)N_x$ and Π the projection matrix onto Ker(A) introduced in the previous part. Let us denote $U(t) = (\tilde{\mathbf{f}}(t), \mathbf{B}(t), \mathbf{E}_1(t), \mathbf{E}_2(t))$, where $\tilde{\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

$$ \begin{aligned} \vec{v}_{\ell} &\in \mathbb{R}^{N_{v_{\ell}}}, \quad (\vec{v}_{\ell})_j = v_{\ell,j} = -v_{\ell,\max} + j \Delta v_{\ell}, \quad j = 1, \dots, N_{v_{\ell}} \text{ and } \ell = 1, 2, \\ \vec{f}_{\star, j_2} &\in \mathbb{R}^{(k+1)N_x N_{v_1}} \quad \text{for } j_2 = 1, \dots, N_{v_2}. \end{aligned} $$
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 (25) can be recast as

$$ \partial_t U = LU + N(U), \qquad (26) $$

with

$$ U = (\vec{\mathbf{f}}_{\star,1}, \vec{\mathbf{f}}_{\star,2}, \dots, \vec{\mathbf{f}}_{\star,N_{v_2}}, \mathbf{B}, \mathbf{E}_2, \mathbf{E}_1)^T, \qquad \vec{\mathbf{f}}_{\star,j_2} \in \mathbb{R}^{\tilde{N}_1}, \forall j_2 = 1, \dots, N_{v_2} \\ N(U) = -(\tilde{\mathcal{F}}\mathcal{D}\vec{\mathbf{f}})_{\star,1}, -(\tilde{\mathcal{F}}\mathcal{D}\vec{\mathbf{f}})_{\star,2}, \dots, -(\tilde{\mathcal{F}}\mathcal{D}\vec{\mathbf{f}})_{\star,N_{v_2}}, \mathbf{0}, \mathbf{0}, \mathbf{0})^T, \qquad (27) $$

$$ L = \begin{pmatrix} \operatorname{diag}(\vec{v}_1) \otimes A & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \cdots \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \\ \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \operatorname{diag}(\vec{v}_1) \otimes A & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \cdots \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \vdots & \vdots & \vdots \\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \cdots & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \operatorname{diag}(\vec{v}_1) \otimes A & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \\ \mathbf{0}_{\tilde{N}, \tilde{N}_1} & \cdots & \mathbf{0}_{\tilde{N}, \tilde{N}_1} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} \\ \mathcal{E}_{2,1} & \mathcal{E}_{2,2} & \cdots & \mathcal{E}_{2,N_{v_2}} & \mathbf{A} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} \\ \mathcal{E}_1 & \mathcal{E}_1 & \cdots & \mathcal{E}_1 & \mathbf{0}_{\tilde{N}, \tilde{N}_1} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} \end{pmatrix} \qquad (28) $$

where $\operatorname{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, $\operatorname{diag}(\vec{v}_1) \otimes A \in \mathbb{M}_{\tilde{N}_1, \tilde{N}_1}$ and $\mathbf{0}_{m,n}$ is the zero matrix with $m$ lines and $n$ columns. Moreover, the matrices $\mathcal{E}_1, \mathcal{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

$$ \begin{aligned} \mathcal{E}_1 &= -\Delta v_1 \Delta v_2 \vec{v}_1 \otimes (\mathbf{1} - \Pi), \\ \mathcal{E}_{2,j_2} &= -\Delta v_1 \Delta v_2 v_{2,j_2} \mathbf{1} \otimes (\mathbf{1} - \Pi) \quad \text{with} \quad \mathbf{1} = (1, \dots, 1) \in \mathbb{R}^{N_{v_1}}. \end{aligned} $$

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 (26) 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 t L) U^n + \Delta t \exp(\Delta t L) N(U^n), \qquad (29) $$

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

**Proposition 4.** The exponential $\exp(L\Delta t)$ of the matrix $L$ given by (28) is given by

$$ \begin{pmatrix} e^{\Delta t \operatorname{diag}(\vec{v}_1) \otimes A} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \cdots \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \\ \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & e^{\Delta t \operatorname{diag}(\vec{v}_1) \otimes A} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \cdots \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \vdots & \vdots & \vdots \\ \vdots & \ddots & \vdots & \vdots & \vdots & \vdots \\ \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \cdots & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & e^{\Delta t \operatorname{diag}(\vec{v}_1) \otimes A} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} & \mathbf{0}_{\tilde{N}_1, \tilde{N}_1} \\ e\mathcal{B}_1 & e\mathcal{B}_2 & \cdots & e\mathcal{B}_{N_{v_2}} & \frac{e^{A\Delta t+e^{-A\Delta t}}}{2} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} \\ e\mathcal{E}_{2,1} & e\mathcal{E}_{2,2} & \cdots & e\mathcal{E}_{2,N_{v_2}} & \frac{e^{A\Delta t-e^{-A\Delta t}}}{2} & \mathbf{0}_{\tilde{N}, \tilde{N}_1} \\ e\mathcal{E}_1 & e\mathcal{E}_1 & \cdots & e\mathcal{E}_1 & \mathbf{0}_{\tilde{N}, \tilde{N}_1} & \mathbf{1} \end{pmatrix} $$
where $e\mathcal{B}_{j_2}, e\mathcal{E}_{2,j_2}, e\mathcal{E}_1 \in \mathbb{M}_{\tilde{N}, \tilde{N}_1}(\mathbb{R})$ for $j_2 = 1, \dots, N_{v_2}$ are given by

$$e\mathcal{B}_{j_2} = \Delta v_1 \Delta v_2 v_{2,j_2} (A + \Pi)^{-1} \vec{\alpha},$$

$$e\mathcal{E}_{2,j_2} = \Delta v_1 \Delta v_2 v_{2,j_2} (A + \Pi)^{-1} \vec{\beta},$$

$$e\mathcal{E}_1 = \Delta v_1 \Delta v_2 (A + \Pi)^{-1} (1 \otimes \mathbf{1} - e^{\Delta t \mathrm{diag}(\vec{v_1})} \otimes A)$$

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

$$\vec{\alpha} = [\alpha_1, \alpha_2, \dots, \alpha_{N_{v_1}}] \in \mathbb{M}_{\tilde{N}, \tilde{N}_1}(\mathbb{R}) \text{ with } \alpha_{j_1} = \left[ \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 t v_{1,j_1}}}{(1 - v_{1,j_1}^2)} \right],$$

$$\vec{\beta} = [\beta_1, \beta_2, \dots, \beta_{N_{v_1}}] \in \mathbb{M}_{\tilde{N}, \tilde{N}_1}(\mathbb{R}) \text{ with } \beta_{j_1} = \left[ \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 t v_{1,j_1}}}{(1 - v_{1,j_1}^2)} \right].$$

*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 \mathrm{diag}(\vec{v_1})} \otimes A$ ($N_{v_2}$ times). Second, the 3 × 3 right bottom block corresponds to the homogeneous Maxwell equations. Its exponential can be computed and is equal to

$$\exp \begin{pmatrix} \mathbf{0}_{\tilde{N}, \tilde{N}} & A \Delta t & \mathbf{0}_{\tilde{N}, \tilde{N}} \\ A \Delta t & \mathbf{0}_{\tilde{N}, \tilde{N}} & \mathbf{0}_{\tilde{N}, \tilde{N}} \\ \mathbf{0}_{\tilde{N}, \tilde{N}} & \mathbf{0}_{\tilde{N}, \tilde{N}} & \mathbf{0}_{\tilde{N}, \tilde{N}} \end{pmatrix} = \begin{pmatrix} \frac{\exp(A \Delta t) + \exp(-A \Delta t)}{2} & \frac{\exp(A \Delta t) - \exp(-A \Delta t)}{2} & \mathbf{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} & \mathbf{0}_{\tilde{N}, \tilde{N}} \\ \mathbf{0}_{\tilde{N}, \tilde{N}} & \mathbf{0}_{\tilde{N}, \tilde{N}} & \mathbf{1} \end{pmatrix}.$$

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

**Computation of $e\mathcal{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})$

$$\mathbf{E}_1(t^{n+1}) = \mathbf{E}_1(t^n) - \sum_{j_1, j_2} \int_{t^n}^{t^{n+1}} (\mathbf{1} - \mathbf{\Pi}) e^{v_{1,j_1} A(t-t^n)} dt v_{1,j_1} \mathbf{f}_{j_1, j_2}(t^n) \Delta v_1 \Delta v_2 \\ = \mathbf{E}_1(t^n) + \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{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\mathcal{B}_{j_2}, e\mathcal{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 us write down the equations for $(\mathbf{E}_2, \mathbf{B})(t)$ with $t \in [t^n, t^{n+1}]$

$$\frac{d}{dt} \mathbf{E}_2(t) = A \mathbf{B}(t) - \sum_{j_1, j_2} (\mathbf{1} - \mathbf{\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$$

$$\frac{d}{dt} \mathbf{B}(t) = 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

$$M = \begin{pmatrix} \mathbf{0}_{\tilde{N}, \tilde{N}} & A \\ A & \mathbf{0}_{\tilde{N}, \tilde{N}} \end{pmatrix}, \quad R(t) = \begin{pmatrix} R_1(t) \\ \mathbf{0}_{\tilde{N}, 1} \end{pmatrix}, \\ R_1(t) = - \sum_{j_1, j_2} (\mathbf{1} - \mathbf{\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 \quad (30)$$
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. \quad (31)$$

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

$$e^{M\Delta t} = \frac{1}{2} \begin{pmatrix} 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{pmatrix}.$$

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

$$\begin{aligned}
\int_{t^n}^{t^{n+1}} e^{-M(t-t^{n+1})} R(t)dt &= \frac{1}{2} \left( \begin{array}{c} \int_{t^n}^{t^{n+1}} \left[ e^{-A(t-t^{n+1})} + e^{A(t-t^{n+1})} \right] R_1(t)dt \\ \int_{t^n}^{t^{n+1}} \left[ e^{-A(t-t^{n+1})} - e^{A(t-t^{n+1})} \right] R_1(t)dt \end{array} \right) \\
&= -\frac{1}{2} \left( \begin{array}{c} \int_{t^n}^{t^{n+1}} \sum_{j_1,j_2} \left[ e^{-A(t-t^{n+1})} + e^{A(t-t^{n+1})} \right] e^{v_{1,j_1}A(t-t^n)} v_{2,j_1} f_{k,j_1,j_2}^n \Delta v_1 \Delta v_2 dt \\ \int_{t^n}^{t^{n+1}} \sum_{j_1,j_2} \left[ e^{-A(t-t^{n+1})} - e^{A(t-t^{n+1})} \right] e^{v_{1,j_1}A(t-t^n)} v_{2,j_1} f_{k,j_1,j_2}^n \Delta v_1 \Delta v_2 dt \end{array} \right) \\
&= - \left( \frac{\sum_{j_1,j_2} \mathcal{I}_1 v_{2,j_1} f_{k,j_1,j_2}^n \Delta v_1 \Delta v_2}{\sum_{j_1,j_2} \mathcal{I}_2 v_{2,j_1} f_{k,j_1,j_2}^n \Delta v_1 \Delta v_2} \right)
\end{aligned}$$

where $\mathcal{I}_1, \mathcal{I}_2$ are given by

$$\begin{aligned}
\mathcal{I}_1 &= \frac{1}{2} \int_{t^n}^{t^{n+1}} \left[ e^{-A(t-t^{n+1})} + e^{A(t-t^{n+1})} \right] e^{v_{1,j_1} A(t-t^n)} dt \\
&= \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 t v_{1,j_1}}}{(1-v_{1,j_1}^2)},
\end{aligned}$$

$$\begin{aligned}
\mathcal{I}_2 &= \frac{1}{2} \int_{t^n}^{t^{n+1}} \left[ e^{-A(t-t^{n+1})} - e^{A(t-t^{n+1})} \right] e^{v_{1,j_1} A(t-t^n)} dt \\
&= \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 t v_{1,j_1}}}{(1-v_{1,j_1}^2)}.
\end{aligned}$$

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

$$\mathbf{E}_2(t^{n+1}) = \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)$$

$$\mathbf{B}(t^{n+1}) = \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

$$\beta_{j_1} = \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 t v_{1,j_1}}}{(1-v_{1,j_1}^2)}$$

$$\alpha_{j_1} = \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 t v_{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 (29) where the exponential is given in Prop 4 (and its generalization to high order Lawson Runge-Kutta) satisfied by $U^n = (\tilde{\mathbf{f}}, \mathbf{B}, \mathbf{E}_2, \mathbf{E}_1)^n$ preserves the following Poisson equation

$$ (A + \Pi)\mathbf{E}_1^n = -\sum_{j_1, j_2} (\mathbf{1} - \Pi)\mathbf{f}_{j_1, j_2}^n \Delta v_1 \Delta v_2, \quad \forall n \in \mathbb{N}^*, $$

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 $\mathbf{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} (\mathbf{1} - \Pi)\mathbf{f}_{j_1, j_2}^0 \Delta v_1 \Delta v_2$ holds at the initial time.

Next, from the scheme (29) 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 t e^{v_{1,j_1} A \Delta t} (\tilde{\mathcal{F}} \mathcal{D} \mathbf{f})_{j_1, j_2}^n. $$

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

$$
\begin{align*}
\mathbf{E}_1^{n+1} &= \mathbf{E}_1^n + \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{1} - e^{v_{1,j_1} A \Delta t}) \mathbf{f}_{j_1, j_2}^n \\
&\quad - \Delta t \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{1} - e^{v_{1,j_1} A \Delta t}) (\tilde{\mathcal{F}} \mathcal{D} \mathbf{f})_{j_1, j_2}^n \\
&= \mathbf{E}_1^n + \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{1} - \Pi) (\mathbf{1} - e^{v_{1,j_1} A \Delta t}) \mathbf{f}_{j_1, j_2}^n \\
&\quad - \Delta t \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{1} - \Pi) (\mathbf{1} - e^{v_{1,j_1} A \Delta t}) (\tilde{\mathcal{F}} \mathcal{D} \mathbf{f})_{j_1, j_2}^n \\
&= \mathbf{E}_1^n + \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{1} - \Pi) (\mathbf{1} - e^{v_{1,j_1} A \Delta t}) \mathbf{f}_{j_1, j_2}^n \\
&\quad - \Delta t \Delta v_1 \Delta v_2 (A + \Pi)^{-1} (\mathbf{1} - \Pi) \sum_{j_1, j_2} (\tilde{\mathcal{F}} \mathcal{D} \mathbf{f})_{j_1, j_2}^n \\
&\quad + \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} (\mathbf{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}) \\
&= \mathbf{E}_1^n + \Delta v_1 \Delta v_2 (A + \Pi)^{-1} \sum_{j_1, j_2} ((\mathbf{1} - \Pi) \mathbf{f}_{j_1, j_2}^n - (\mathbf{1} - \Pi) \mathbf{f}_{j_1, j_2}^{n+1}).
\end{align*}
$$

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. $\square$

# 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 \ge 0$, $(x, y) \in [0, L_x] \times [0, L_y]$ and $(v_1, v_2) \in$
$$
\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 \, dv_1 \, dv_2 + \bar{J}_1, \quad \partial_t E_2 = -\partial_x B - \int_{\mathbb{R}^2} v_2 f \, dv_1 \, dv_2 + \bar{J}_2, \\
& \partial_x E_1 + \partial_y E_2 = \int_{\mathbb{R}^2} f \, dv_1 \, dv_2 - \bar{\rho}, \quad \partial_x B + \partial_y B = 0,
\end{aligned}
\right.
\qquad (32)
$$

with initial conditions ($f_0(x, y, v_1, v_2)$, $E_1^0(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 dv_1 dv_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 f dx dy dv_1 dv_2$, $\bar{\rho} = \frac{1}{L_x L_y} \int_{L_x \times L_y} \int_{\mathbb{R}^2} f dx dy dv_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, \dots, N_x$ and $j = 1, 2, \dots, N_y$):

$$
\begin{align}
& \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] \nonumber \\
& - \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] \nonumber \\
& + \sum_{m=0}^{k} \sum_{n=0}^{k} \left[ (v_{j_1} [\{f_h(t, x, y, v_{j_1}, v_{j_2})\} \xi^\ell]_{i-\frac{1}{2}}^{i+\frac{1}{2}}, \xi^s)_{I_j} + (v_{j_2} [\{f_h(t, x, y, v_{j_1}, v_{j_2})\} \xi^s]_{j-\frac{1}{2}}^{j+\frac{1}{2}}, \xi^\ell)_{I_i} \right] = 0, \tag{33}
\end{align}
$$

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 $f_{j_1, j_2} \in \mathbb{R}^{(k+1)^2 N_x N_y}$ the vector containing the degree of freedom $f_{i, j}^{m, n}$

$$\mathbf{f}_{j_1, j_2} = [f_{1,1}^{0,0}, f_{1,1}^{1,0}, \dots, f_{1,1}^{k,0}, \dots, f_{N_x,1}^{k,0}, f_{1,1}^{0,1}, \dots, f_{N_x, N_y}^{k,k}]_{j_1, j_2}^T \quad (34)$$

and $\mathbf{E}_1, \mathbf{E}_2, \mathbf{B} \in \mathbb{R}^{(k+1)^2 N_x N_y}$ the vectors (defined as (34)) containing the DG degree of freedom of the electromagnetic fields $(E_1, E_2, B)$. Finally, we introduce $(\mathcal{F}\mathcal{D}\mathbf{f})_{j_1, j_2}$ the DG approximation of the nonlinear term $(\mathcal{F} \cdot \nabla_v f)(v_{j_1}, v_{j_2})$, with $\mathcal{F} = (E_1 + Bv_2, E_2 - Bv_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} &= (\mathbf{1}_y \otimes v_{1,j_1} A^x) \mathbf{f}_{j_1, j_2} + (v_{2,j_2} A^y \otimes \mathbf{1}_x) \mathbf{f}_{j_1, j_2} - (\tilde{\mathcal{F}} \mathcal{D} \mathbf{f})_{j_1, j_2} \\
\partial_t \mathbf{B} &= -(A^y \otimes \mathbf{1}_x) \mathbf{E}_1 + (\mathbf{1}_y \otimes A^x) \mathbf{E}_2, \\
\partial_t \mathbf{E}_1 &= -(A^y \otimes \mathbf{1}_x) \mathbf{B} - \sum_{j_1, j_2} v_{1,j_1} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2, \\
\partial_t \mathbf{E}_2 &= (\mathbf{1}_y \otimes A^x) \mathbf{B} - \sum_{j_1, j_2} v_{2,j_2} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2,
\end{aligned}
\right.
\quad (35)
$$

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, $\mathbf{P} = (\mathbf{1}_y - \Pi_y) \otimes (\mathbf{1}_x - \Pi_x)$ with $\mathbf{1}_x$ (resp. $\mathbf{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 $\text{Ker}(A^x)$ (resp. $\text{Ker}(A^y)$).

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

**Proposition 6.** The semi-discretized system (35) satisfied by $(\mathbf{f}_{j_1,j_2}, \mathbf{B}, \mathbf{E}_1, \mathbf{E}_2)(t)$ preserves the following discretized Poisson equation (provided that it is satisfied at time $t=0$)

$$ (\mathbf{1}_y \otimes A^x)\mathbf{E}_1(t) + (A^y \otimes \mathbf{1}_x)\mathbf{E}_2(t) = - \sum_{j_1, j_2} \mathbf{P} \mathbf{f}_{j_1, j_2}(t) \Delta v_1 \Delta v_2, \quad \forall t > 0. $$

*Proof.* Let us derive with respect to time the left hand side of the equality to get

$$
\begin{align*}
\partial_t((\mathbf{1}_y \otimes A^x)\mathbf{E}_1 + (A^y \otimes \mathbf{1}_x)\mathbf{E}_2) &= (\mathbf{1}_y \otimes A^x)\partial_t\mathbf{E}_1 + (A^y \otimes \mathbf{1}_x)\partial_t\mathbf{E}_2 \\
&= (\mathbf{1}_y \otimes A^x)(-(A^y \otimes \mathbf{1}_x)\mathbf{B} - \sum_{j_1, j_2} v_{1,j_1} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2) \\
&\quad + (A^y \otimes \mathbf{1}_x)((\mathbf{1}_y \otimes A^x)\mathbf{B} - \sum_{j_1, j_2} v_{2,j_2} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2) \\
&= [(A^y \otimes \mathbf{1}_x)(\mathbf{1}_y \otimes A^x) - (\mathbf{1}_y \otimes A^x)(A^y \otimes \mathbf{1}_x)]\mathbf{B} \\
&\quad - (\mathbf{1}_y \otimes A^x) \sum_{j_1, j_2} v_{1,j_1} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2) - (A^y \otimes \mathbf{1}_x) \sum_{j_1, j_2} v_{2,j_2} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2) \\
&= -\sum_{j_1, j_2} \left(\partial_t(\mathbf{P} \mathbf{f}_{j_1, j_2}) + \mathbf{P}(\tilde{\mathcal{F}} \mathcal{D} \mathbf{f})_{j_1, j_2}\right) \Delta v_1 \Delta v_2 = -\partial_t\left(\sum_{j_1, j_2} \mathbf{P} \mathbf{f}_{j_1, j_2} \Delta v_1 \Delta v_2\right),
\end{align*}
$$

where we used the identities

$$ (A^y \otimes \mathbf{1}_x)(\mathbf{1}_y \otimes A^x) = (A^y \mathbf{1}_y) \otimes (\mathbf{1}_x A^x) = A^y \otimes A^x $$

$$ (\mathbf{1}_y \otimes A^x)(A^y \otimes \mathbf{1}_x) = (\mathbf{1}_y A^y) \otimes (A^x \mathbf{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. $\square$

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 $\mathbf{f}_{j_1,j_2}$ linear part (which corresponds to the $(x, y)$ transport). Indeed, we observe from (35) that this linear part writes

$$ \partial_t \mathbf{f}_{j_1, j_2} = [(\mathbf{1}_y \otimes v_{1,j_1} A^x) + (v_{2,j_2} A^y \otimes \mathbf{1}_x)] \mathbf{f}_{j_1, j_2} = [v_{2,j_2} A^y \otimes v_{1,j_1} A^x] \mathbf{f}_{j_1, j_2}, $$
where we used the definition of the Kronecker sum $\oplus$. The exact solution can be then written as

$$ \mathbf{f}_{j_1, j_2}(t) = \exp \left( (v_{2,j_2} A^y \oplus v_{1,j_1} A^x) t \right) \mathbf{f}_{j_1, j_2}(0) = \exp \left( v_{2,j_2} t A^y \oplus v_{1,j_1} t A^x \right) \mathbf{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

$$ \mathbf{f}_{j_1, j_2}(t) = \exp(v_{2,j_2} t A^y) \otimes \exp(v_{1,j_1} t A^x) \mathbf{f}_{j_1, j_2}(0), $$

which can be recast using the vec operation as

$$ \mathbf{t}_{j_1, j_2}(t) = \exp(v_{2,j_2} t A^y) \mathbf{t}_{j_1, j_2}(0) \exp(\Delta v_{1,j_1} t (A^x)^T), \quad (36) $$

where $\mathbf{f}_{j_1, j_2}(s) = \text{vec}(\mathbf{t}_{j_1, j_2}(s))$ (for $s = 0, t$ denotes the vectorization operation which takes the matrix $\mathbf{t}_{j_1, j_2} \in \mathbb{R}^{(k+1)N_x (k+1)N_y}$ as entry and gives the vector $\mathbf{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 [7, 17, 29].

**Remark 5.** The semi-discretized Vlasov-Maxwell system (35) 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 $\mathbf{f}$ and then requires the calculation of $\exp(v_{2,j_2} t A^y \oplus v_{1,j_1} t A^x)$ which can be performed efficiently thanks to (36). The update of the electric field is performed using the Poisson equation thanks to the updated $\mathbf{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 (32) and we will see that in this case, it will be possible to compute explicitly 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 $\mathcal{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} i k_x + v_{2,j_2} i k_y) \hat{f}_{k_x, k_y, j_1, j_2} + (\widehat{\mathcal{F}Df})_{k_x, k_y, j_1, j_2} &= 0, \\ \partial_t \hat{B}_{k_x, k_y} &= i k_y \hat{E}_{1, k_x, k_y} - i k_x \hat{E}_{2, k_x, k_y}, \\ \partial_t \hat{E}_{1, k_x, k_y} &= i k_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} &= -i k_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. \quad (37) $$

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) + iky \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 (37), we have a similar proposition as Prop (6) in this Fourier case.
**Proposition 7.** The semi-discretized system (37) 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 (provided that it is satisfied at time $t=0$)

$$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, \quad (k_x, k_y) \neq (0,0), \forall 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

$$\begin{aligned} \partial_t (ik_x \hat{E}_{1,k_x,k_y} + ik_y \hat{E}_{2,k_x,k_y}) &= ik_x (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) \\ &\quad + ik_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) \\ &= \sum_{j_1,j_2} \Delta v_1 \Delta v_2 (\partial_t \hat{f}_{k_x,k_y,j_1,j_2} + (\widehat{\mathcal{FD}}f)_{k_x,k_y,j_1,j_2}) = \partial_t (\sum_{j_1,j_2} \hat{f}_{k_x,k_y,j_1,j_2} \Delta v_1 \Delta v_2), \end{aligned}$$

where we used the summation on $j_1, j_2$ of the Vlasov equation together with the fact that $\mathcal{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. $\square$

To derive a fully discrete scheme, we introduce as previously the vector $\hat{f}_{k_x,k_y} \in \mathbb{C}^{N_{v_1} N_{v_2}}$ and denote $U(t) = (\hat{f}, \hat{B}, \hat{E}_2, \hat{E}_1)_{k_x,k_y}(t) \in \mathbb{C}^{N_{v_1} N_{v_2}+3}$, then the system (37) can be rewritten as

$$\partial_t U = LU + N(U), \quad (38)$$

where $L$ is given by

$$\left( \begin{array}{c} -i(k_x \mathrm{diag}(\vec{v}_1) + k_y v_{2,1}) \quad \mathbf{0}_{N_{v_1}, N_{v_1}} \quad \dots \quad \mathbf{0}_{N_{v_1}, N_{v_1}} \quad \mathbf{0}_{N_{v_1}, 1} \quad \dots \quad \mathbf{0}_{N_{v_1}, 1} \\ \mathbf{0}_{N_{v_1}, N_{v_1}} \quad -i(k_x \mathrm{diag}(\vec{v}_1) + k_y v_{2,2}) \quad \dots \quad \mathbf{0}_{N_{v_1}, N_{v_1}} \quad \vdots \quad \dots \quad \vdots \\ \vdots \quad \dots \quad \dots \quad \mathbf{0}_{N_{v_1}, N_{v_1}} \quad \dots \quad \dots \quad \dots \\ \mathbf{0}_{N_{v_1}, N_{v_1}} \quad \dots \quad \mathbf{0}_{N_{v_1}, N_{v_1}} \quad -i(k_x \mathrm{diag}(\vec{v}_1) + k_y v_{2,N_{v_2}}) \quad \mathbf{0}_{N_{v_1}, 1} \quad \dots \quad \mathbf{0}_{N_{v_1}, 1} \\ \mathbf{0}_{1, N_{v_1}} \quad \dots \quad \mathbf{0}_{1, N_{v_1}} \quad \mathbf{0}_{1, N_{v_1}} \quad \mathbf{0}_{1, N_{v_1}} \quad 0 \quad -ik_x \quad ik_y \\ -\Delta v_1 \Delta v_2 v_{2,1} \mathbf{1} \quad -\Delta v_1 \Delta v_2 v_{2,2} \mathbf{1} \quad \dots \quad -\Delta v_1 \Delta v_2 v_{2,N_{v_2}} \mathbf{1} \quad -ik_x \quad 0 \quad 0 \\ -\Delta v_1 \Delta v_2 \vec{v}_1 \quad -\Delta v_1 \Delta v_2 \vec{v}_1 \quad \dots \quad -\Delta v_1 \Delta v_2 \vec{v}_1 \quad ik_y \quad 0 \quad 0 \end{array} \right) \quad (39)$$

and

$$\hat{\mathbf{f}}_{k_x,k_y} = (\hat{f}_{k_x,k_y,1,1}, \dots, \hat{f}_{k_x,k_y,N_{v_1},1}, \dots, \hat{f}_{k_x,k_y,N_{v_1},N_{v_2}})^T \in \mathbb{C}^{N_{v_1} N_{v_2}},$$

$$N(U) = ((\widehat{\mathcal{FD}}f)_{k_x,k_y}, 0, 0, 0)^T \in \mathbb{C}^{N_{v_1} N_{v_2}+3},$$

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

We now study the time discretization of (38) 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), \quad (40)$$

that approximates the ODE (38) preserves a discrete Poisson equation. This is the object of the following proposition.
**Proposition 8.** The Lawson scheme (40) satisfied by $U^n = (\hat{f}, \hat{B}, \hat{E}_2, \hat{E}_1)_{k_x, k_y}^n$ preserves the following Poisson equation (provided that it is satisfied at the initial time $n = 0$)

$$ 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}_{k_x,k_y,j_1,j_2}^n \Delta v_1 \Delta v_2, \quad (k_x, k_y) \neq (0,0), \forall n \in \mathbb{N}^*, $$

*Proof.* First, we need to know the shape of $\exp(\Delta t L)$. From the one-dimensional calculations and from [5], we have the representation of $e^{L\Delta t}$

$$ \left( \begin{array}{ccccccc} e^{-i\Delta t(k_x v_1 + k_y v_{2,1})} & \mathbf{0}_{N_{v_1},N_{v_1}} & \mathbf{0}_{N_{v_1},N_{v_1}} & \cdots & \mathbf{0}_{N_{v_1},N_{v_1}} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} \\ \mathbf{0}_{N_{v_1},N_{v_1}} & e^{-i\Delta t(k_x v_1 + k_y v_{2,2})} & \mathbf{0}_{N_{v_1},N_{v_1}} & \cdots & \mathbf{0}_{N_{v_1},N_{v_1}} & \vdots & \vdots & \vdots \\ \vdots & \cdots & \vdots & \cdots & \mathbf{0}_{N_{v_1},N_{v_1}} & \vdots & \vdots & \vdots \\ \mathbf{0}_{N_{v_1},N_{v_1}} & \cdots & \mathbf{0}_{N_{v_1},N_{v_1}} & \cdots & e^{-i\Delta t(k_x v_1 + k_y v_{2,N_{v_2}})} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} \\ e\mathcal{B}_{k_x,k_y,\star,1} & e\mathcal{B}_{k_x,k_y,\star,2} & \cdots & \cdots & e\mathcal{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\mathcal{E}_{2,k_x,k_y,\star,1} & e\mathcal{E}_{2,k_x,k_y,\star,2} & \cdots & \cdots & e\mathcal{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\mathcal{E}_{1,k_x,k_y,\star,1} & e\mathcal{E}_{1,k_x,k_y,\star,2} & \cdots & \cdots & e\mathcal{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 \begin{pmatrix} 0 & -ik_x & ik_y \\ -ik_x & 0 & 0 \\ ik_y & 0 & 0 \end{pmatrix} \right) = \begin{pmatrix} \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_y^2}{|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{pmatrix}, $$

and the vectors $(e\mathcal{B}, e\mathcal{E}_2, e\mathcal{E}_1)_{k_x,k_y,\star,j_2} \in 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{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 \quad (41) $$

which can be inserted in the Maxwell part to compute the vectors $(e\mathcal{B}, e\mathcal{E}_2, e\mathcal{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 = \mathcal{J}V + R \text{ with } R(t) = \begin{pmatrix} 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{pmatrix} \text{ and } \mathcal{J} = \begin{pmatrix} 0 & -ik_x & ik_y \\ -ik_x & 0 & 0 \\ ik_y & 0 & 0 \end{pmatrix}. $$

The variation of constant formula gives

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

The calculations for the integral term involve the following integral term

$$ \begin{aligned} \mathcal{I}_1(\hat{f}^n) &= \int_{t^n}^{t^{n+1}} \left( \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) \right) dt \\ &= \sum_{j_1,j_2} \left( \frac{(-ik_x v_{2,j_2} + ik_y v_{1,j_1})}{|k|} C \right) \hat{f}_{k_x,k_y,j_1,j_2}^n \Delta v_1 \Delta v_2, \\ &= \sum_{j_1,j_2} e \mathcal{B}_{k_x,k_y,j_1,j_2} \hat{f}_{k_x,k_y,j_1,j_2}^n, \end{aligned} $$
$$
\begin{align*}
\mathcal{I}_2(\hat{\mathbf{f}}^n) &= \int_{t^n}^{t^{n+1}} \left( \frac{k_x^2 \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) \right) dt \\
&= \sum_{j_1, j_2} \left( v_{2,j_2} \frac{k_x^2 \mathcal{A} + k_y^2 \mathcal{B}}{|k|^2} + v_{1,j_1} \frac{k_x k_y (\mathcal{B} - \mathcal{A})}{|k|^2} \right) \hat{f}_{k_x, k_y, j_1, j_2}^n \Delta v_1 \Delta v_2, \\
&= \sum_{j_1, j_2} e^{\mathcal{E}_2(k_x, k_y, j_1, j_2)} \hat{f}_{k_x, k_y, j_1, j_2}^n,
\end{align*}
$$

$$
\begin{align*}
\mathcal{I}_3(\hat{\mathbf{f}}^n) &= \int_{t^n}^{t^{n+1}} \left( \frac{k_x k_y (1 - \cos(|k|(t^{n+1} - t)))}{|k|^2} R_2(t) + \frac{k_y^2 \cos(|k|(t^{n+1} - t)) + k_x^2}{|k|^2} R_3(t) \right) dt \\
&= \sum_{j_1, j_2} \left( v_{2,j_2} \frac{k_x k_y (\mathcal{B} - \mathcal{A})}{|k|^2} + v_{1,j_1} \frac{k_y^2 \mathcal{A} + k_x^2 \mathcal{B}}{|k|^2} \right) \hat{f}_{k_x, k_y, j_1, j_2}^n \Delta v_1 \Delta v_2, \\
&= \sum_{j_1, j_2} e^{\mathcal{E}_1(k_x, k_y, j_1, j_2)} \hat{f}_{k_x, k_y, j_1, j_2}^n,
\end{align*}
$$

where the time integrals are

$$
\mathcal{A} = \int_{t^n}^{t^{n+1}} \cos(|k|(t^{n+1} - t))e^{-i(t-t^n)(k \cdot v_j)} dt,
$$

$$
\mathcal{B} = \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 t k \cdot v_j}),
$$

$$
\mathcal{C} = \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

$$
\hat{E}_2^{n+1} = (e^{\Delta t \mathcal{J}} V^n)_2 + \mathcal{I}_2(\hat{\mathbf{f}}^n - \Delta t (\widehat{\mathcal{F}\mathcal{D}\mathbf{f}^n})), \quad \hat{E}_1^{n+1} = (e^{\Delta t \mathcal{J}} V^n)_3 + \mathcal{I}_3(\hat{\mathbf{f}}^n - \Delta t (\widehat{\mathcal{F}\mathcal{D}\hat{\mathbf{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

$$
\begin{align*}
ik_y (e^{\Delta t \mathcal{J}} V^n)_2 + ik_x (e^{\Delta t \mathcal{J}} V^n)_3 &= ik_y \left( \frac{-ik_x \sin(|k|\Delta t)}{|k|} \hat{B}^n + \frac{k_x^2 \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 \right) \\
&\quad + ik_x \left( \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_y^2 \cos(|k|\Delta t) + k_x^2}{|k|^2} \hat{E}_1^n \right) = ik_y \hat{E}_2^n + ik_x \hat{E}_1^n.
\end{align*}
$$
Then, we have

$$
\begin{align*}
ik_y \hat{E}_2^{n+1} + ik_x \hat{E}_1^{n+1} &= ik_y \hat{E}_2^n + ik_x \hat{E}_1^n + ik_y \mathcal{I}_2(\hat{f}^n - \Delta t(\widehat{\mathcal{D}f^n})) + ik_x \mathcal{I}_3(\hat{f}^n - \Delta t(\widehat{\mathcal{D}f^n})) \\
&= 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{\mathcal{D}f^n}))\mathcal{B}\Delta v_2\Delta v_1 \\
&\quad + ik_x \sum_{j_1,j_2} v_{1,j_1}(\hat{f}^n - \Delta t(\widehat{\mathcal{D}f^n}))\mathcal{B}\Delta v_2\Delta v_1 \\
&= 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{\mathcal{D}f^n})) \frac{1}{ik \cdot v_j} (1 - e^{-i\Delta t \cdot v_j}) \\
&= ik_y \hat{E}_2^n + ik_x \hat{E}_1^n + \Delta v_2\Delta v_1 \left[ \sum_{j_1,j_2} (1 - e^{-i\Delta t \cdot v_j})\hat{f}^n - \Delta t \sum_{j_1,j_2} (\widehat{\mathcal{D}f^n}) + \Delta t \sum_{j_1,j_2} e^{-i\Delta t \cdot v_j} (\widehat{\mathcal{D}f^n}) \right] \\
&= ik_y \hat{E}_2^n + ik_x \hat{E}_1^n + \Delta v_2\Delta v_1 \left[ \sum_{j_1,j_2} (1 - e^{-i\Delta t \cdot v_j})\hat{f}^n + \sum_{j_1,j_2} e^{-i\Delta t \cdot v_j} \hat{f}^n - \hat{f}^{n+1} \right] \\
&= 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}),
\end{align*}
$$

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

# 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, \quad (x, v) \in [0, 2\pi]^2 \qquad (42) $$

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 an 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 $\mathcal{D}$ is chosen as a 4th order centered finite difference operator CD4), we expect the following estimate

$$ \|\mathbf{u}(t^n) - \mathbf{u}_h^n\|_{L^2} \le 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 [28]. 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 $\mathcal{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.

<table><caption>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$)).</caption><thead><tr><th rowspan="2"></th><th rowspan="2">N<sub>x</sub></th><th colspan="4">central flux</th><th colspan="4">upwind flux</th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>2</sup>-error</th><th>order</th><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>2</sup>-error</th><th>order</th></tr></thead><tbody><tr><td rowspan="5">P<sup>1</sup></td><td>10</td><td>1.62E-01</td><td>-</td><td>3.76E-01</td><td>-</td><td>4.78E-02</td><td>-</td><td>9.51E-02</td><td>-</td></tr><tr><td>20</td><td>7.66E-02</td><td>1.08</td><td>1.86E-01</td><td>1.02</td><td>1.27E-02</td><td>1.91</td><td>2.42E-02</td><td>1.98</td></tr><tr><td>40</td><td>3.70E-02</td><td>1.05</td><td>9.24E-02</td><td>1.01</td><td>3.25E-03</td><td>1.97</td><td>6.07E-03</td><td>1.99</td></tr><tr><td>80</td><td>1.82E-02</td><td>1.03</td><td>4.60E-02</td><td>1.00</td><td>8.23E-04</td><td>1.99</td><td>1.52E-03</td><td>2.00</td></tr><tr><td>160</td><td>9.01E-03</td><td>1.01</td><td>2.30E-02</td><td>1.00</td><td>2.07E-04</td><td>1.99</td><td>3.79E-04</td><td>2.00</td></tr><tr><td rowspan="5">P<sup>2</sup></td><td>10</td><td>2.29E-03</td><td>-</td><td>3.76E-03</td><td>-</td><td>2.52E-03</td><td>-</td><td>4.67E-03</td><td>-</td></tr><tr><td>20</td><td>2.65E-04</td><td>3.11</td><td>4.53E-04</td><td>3.05</td><td>3.07E-04</td><td>3.04</td><td>5.83E-04</td><td>3.00</td></tr><tr><td>40</td><td>3.25E-05</td><td>3.02</td><td>5.63E-05</td><td>3.01</td><td>3.84E-05</td><td>3.00</td><td>7.29E-05</td><td>3.00</td></tr><tr><td>80</td><td>4.05E-06</td><td>3.01</td><td>7.03E-06</td><td>3.00</td><td>4.79E-06</td><td>3.00</td><td>9.11E-06</td><td>3.00</td></tr><tr><td>160</td><td>5.12E-07</td><td>2.98</td><td>8.99E-07</td><td>2.97</td><td>5.95E-07</td><td>3.00</td><td>1.16E-06</td><td>2.98</td></tr></tbody></table>

<table><caption>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$)).</caption><thead><tr><th rowspan="2"></th><th rowspan="2">N<sub>v</sub></th><th colspan="4">central flux</th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>2</sup>-error</th><th>order</th></tr></thead><tbody><tr><td rowspan="5">CD4</td><td>8</td><td>1.18E-02</td><td>-</td><td>5.24E-02</td><td>-</td></tr><tr><td>16</td><td>7.78E-04</td><td>3.92</td><td>3.50E-03</td><td>3.92</td></tr><tr><td>32</td><td>4.93E-05</td><td>3.98</td><td>2.19E-04</td><td>3.98</td></tr><tr><td>64</td><td>3.09E-06</td><td>4.00</td><td>1.37E-05</td><td>4.00</td></tr><tr><td>128</td><td>1.98E-07</td><td>3.97</td><td>8.78E-07</td><td>3.97</td></tr></tbody></table>

## 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 $H(t) - H(0)$ with $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><img src="image_3.png" /><figcaption>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$).</figcaption></figure>

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 particular, 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.

<table><caption>Table 3: Linear transport equation: $L^∞$ 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$)).</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Δt</th><th colspan="4">central flux</th></tr><tr><th>L<sup>∞</sup>-error</th><th>order</th><th>L<sup>2</sup>-error</th><th>order</th></tr></thead><tbody><tr><td rowspan="5">RK(3,3)</td><td>0.1000</td><td>4.17E-05</td><td>-</td><td>1.85E-04</td><td>-</td></tr><tr><td>0.0500</td><td>5.21E-06</td><td>2.99</td><td>2.31E-05</td><td>2.99</td></tr><tr><td>0.0250</td><td>6.51E-07</td><td>3.00</td><td>2.89E-06</td><td>3.00</td></tr><tr><td>0.0125</td><td>8.14E-08</td><td>3.00</td><td>3.61E-07</td><td>3.00</td></tr><tr><td>0.00625</td><td>1.02E-08</td><td>3.00</td><td>4.52E-08</td><td>3.00</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>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$).</figcaption></figure>

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

We consider the Weibel instability [37] by considering 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)), \quad x \in [0, 2\pi/k),$$ $$B(t=0x) = \beta \cos(kx), \quad 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 [37]). 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_{v1} = N_{v2} = 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_{v1} = N_{v2} = 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 $|\mathcal{H}(t) - \mathcal{H}(0)|/\mathcal{H}(0)$ in Figure 4 for Lawson-RK(3,3)-DG CD4 scheme without (left) and with (right) the energy correction step presented in [5] 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><img src="image_5.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

The second test for Vlasov-Maxwell equation we considered is the streaming Weibel instability [8, 11] 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), \quad 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 [8]. 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><img src="image_7.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Vlasov-Maxwell equation (streaming Weibel instability): time evolution of the relative total energy of Lawson-RK(3,3)-DG CD4 without (left) and with (right) energy correction step.</figcaption></figure>

# 6 Conclusion

In this work, we constructed and implemented a new exponential DG method for Vlasov equation, extending the previous works [18, 16, 5] 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.

## 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_j^\ell$ and $u_h(t, x) = \sum_{m=0}^k u_j^m(t) \xi_j^m(x)$ in (2) to get

$$ \sum_{m=0}^{k} \left[ (\partial_t u_j^m(t) \xi^m, \xi^\ell)_{I_j} - (u_j^m(t) \xi^m, \partial_x \xi^\ell)_{I_j} \right] + [\{u_h\}\xi^\ell]_{x_{j-1/2}}^{x_{j+1/2}} = 0, \quad (43) $$

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

$$ \begin{aligned} [\{u_h\}\xi^\ell]_{x_{j-1/2}}^{x_{j+1/2}} &= \frac{1}{2} \sum_{m=0}^{k} \left[ (u_j^m(t)\xi^{m,-}(x_{j+1/2}) + u_{j+1}^m(t)\xi^{m,+}(x_{j+1/2}))\xi^{\ell,-}(x_{j+1/2}) \right. \\ & \qquad \left. - (u_{j-1}^m(t)\xi^{m,-}(x_{j-1/2}) + u_j^m(t)\xi^{m,+}(x_{j-1/2}))\xi^{\ell,+}(x_{j-1/2}) \right] \\ &= \frac{1}{2} \sum_{m=0}^{k} \left[ (u_j^m(t)(1/2)^m + u_{j+1}^m(t)(-1/2)^m)(1/2)^\ell \right. \\ & \qquad \left. - (u_{j-1}^m(t)(1/2)^m + u_j^m(t)(-1/2)^m)(-1/2)^\ell \right]. \end{aligned} $$

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

$$ \Delta x \begin{pmatrix} M \\ M \\ \vdots \\ M \end{pmatrix} \partial_t \begin{pmatrix} \mathbf{u}_1 \\ \mathbf{u}_2 \\ \vdots \\ \mathbf{u}_N \end{pmatrix} - \begin{pmatrix} D_1 & D_2 & \dots & D_3 \\ D_3 & D_1 & D_2 & \vdots \\ \vdots & \vdots & \ddots & D_2 \\ D_2 & \dots & D_3 & D_1 \end{pmatrix} \begin{pmatrix} \mathbf{u}_1 \\ \mathbf{u}_2 \\ \vdots \\ \mathbf{u}_N \end{pmatrix} = \mathbf{0}, \quad (44) $$

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

$$ \begin{aligned} M_{\ell,m} &= \frac{(1/2)^{m+\ell-1}}{m+\ell-1} [1 - (-1)^{m+\ell-1}], \quad (D_2)_{\ell,m} = (-1)^m (1/2)^{m+\ell-1}, \\ (D_1)_{\ell,m} &= (1/2)^{m+\ell-2} \left( \frac{\ell-1}{m+\ell-2} - \frac{1}{2} \right) [1 - (-1)^{m+\ell-2}], \quad (D_3)_{\ell,m} = (-1)^{\ell-1} (1/2)^{m+\ell-1}, \end{aligned} $$

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

$$ \Delta x \mathcal{M} \partial_t \mathbf{u} - \mathcal{D} \mathbf{u} = 0, $$
where $\mathcal{M}$ is a block diagonal mass matrix of size $N(k+1)$, $\mathcal{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), \quad A = (\Delta x \mathcal{M})^{-1} \mathcal{D}, \qquad (45) $$

where $A$ is a block circulant matrix

$$ A = \frac{1}{\Delta x} \begin{pmatrix} M^{-1}D_1 & M^{-1}D_2 & \mathbf{0} & \dots & M^{-1}D_3 \\ M^{-1}D_3 & M^{-1}D_1 & M^{-1}D_2 & \mathbf{0} & \dots \\ \mathbf{0} & \ddots & \ddots & \ddots & M^{-1}D_2 \\ M^{-1}D_2 & \mathbf{0} & \dots & M^{-1}D_3 & M^{-1}D_1 \end{pmatrix} = \frac{1}{\Delta x} A_1, $$

or with the circblock notation $A_1 = \text{circblock}(M^{-1}D_1, M^{-1}D_2, 0, 0, \dots, M^{-1}D_3)$. Since $\mathcal{D}$ and $\mathcal{M}$ are independent of $\Delta x$, so is $A_1$ defined by $A = \Delta x^{-1}A_1$.

## 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,

$$ \|v\| = (\sum_j \|v\|_j^2)^{\frac{1}{2}}, \quad \|v\|_{\infty} = \max_j \|v\|_{\infty,j}, \quad \|v\|_{\Gamma}^2 = \sum_j (|v_{j-1/2}^+|^2 + |v_{j-1/2}^-|^2), $$

where we express the value of $v$ on the left and right limits of the grid point $x_{j+1/2}$ with $v_{j+1/2}^{-}$ and $v_{j+1/2}^{+}$ respectively. Define the jump and the mean of $v$ at $x_{j-1/2}$ as $[v]_{j-1/2} = (v_{j-1/2}^{+} - v_{j-1/2}^{-})$ and $\{v\}_{j-1/2} = (v_{j-1/2}^{+} + v_{j-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 \le C \|v\|^2, \qquad (46) $$

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:

$$ (P_k u, v_h) = (u, v_h), \forall v_h \in V_h \text{ with } v_h(\chi_j(\cdot, t)) \in P^k([-1, 1]). $$

The following lemma states the error of these projections [13].
**Lemma B.2.** Let $P_h^\perp q = q - P_k q$ be 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} \leq ch^{k+1} |q|_{k+1, D}, \quad (47) $$

$$ \|P_h^\perp q\|_{\Gamma} \leq ch^{k+\frac{1}{2}} \|\partial_x^{k+1} q\|, \quad (48) $$

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

$$ \begin{aligned} \frac{1}{2} \frac{d}{dt} \int_{I_j} u_h^2 dx &= -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, \\ &= -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}}}. \\ &= -\frac{a}{2}u_h^+ u_h^-|_{x_{j+\frac{1}{2}}} + \frac{a}{2}u_h^- u_h^+|_{x_{j-\frac{1}{2}}}. \end{aligned} \quad (49) $$

Sum over $j$ of above equation, the $L^2$-stability follows. $\square$

**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} \leq Ch^k, \quad (50) $$

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

$$ \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. \quad (51) $$

Define $e_u = u_h - u = (u_h - P_k u) - (u - P_k u) = \tilde{e}_u - P_k^\perp u$. Then taking $\psi = \tilde{e}_u$, we have

$$ \begin{aligned} & \int_{I_j} \partial_t (\tilde{e}_u) \tilde{e}_u dx \\ &= \int_{I_j} \partial_t (P_k^\perp u) \tilde{e}_u dx - a\{\tilde{e}_u\}|_{x_{j+\frac{1}{2}}} \tilde{e}_u^-|_{x_{j+\frac{1}{2}}} + a\{\tilde{e}_u\}|_{x_{j-\frac{1}{2}}} \tilde{e}_u^+|_{x_{j-\frac{1}{2}}} + \int_{I_j} a \tilde{e}_u (\tilde{e}_u)_x dx \quad (52) \\ & \qquad + a\{P_k^\perp u\}|_{x_{j+\frac{1}{2}}} \tilde{e}_u^-|_{x_{j+\frac{1}{2}}} - a\{P_k^\perp u\}|_{x_{j-\frac{1}{2}}} \tilde{e}_u^+|_{x_{j-\frac{1}{2}}} - \int_{I_j} a P_k^\perp u (\tilde{e}_u)_x dx. \end{aligned} $$
By the definition of the projections and some calculations, the right terms RHS of the error equation (52) become

$$RHS = -\frac{a}{2}\tilde{e}_u^+\tilde{e}_u^-|_{x_{j+\frac{1}{2}}} + \frac{a}{2}\tilde{e}_u^-\tilde{e}_u^+|_{x_{j-\frac{1}{2}}} + a\{P_k^\perp u\}|_{x_{j+\frac{1}{2}}}\tilde{e}_u^-|_{x_{j+\frac{1}{2}}} - a\{P_k^\perp u\}|_{x_{j-\frac{1}{2}}}\tilde{e}_u^+|_{x_{j-\frac{1}{2}}}. \quad (53)$$

Sum over $j$,

$$\frac{1}{2}\frac{d}{dt}\|\tilde{e}_u\|^2 = \sum_j \text{RHS} = -\sum_j a\{P_k^\perp u\}[\tilde{e}_u]|_{x_{j+\frac{1}{2}}}.$$

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

$$\begin{aligned} \frac{1}{2}\frac{d}{dt}\|\tilde{e}_u\|^2 &\le \sum_j \left|a\{P_k^\perp u\}[\tilde{e}_u]|_{x_{j+\frac{1}{2}}}\right| \\ &\le c\|P_k^\perp u\|_{\Gamma}\|\tilde{e}_u\|_{\Gamma} \le ch^{k+\frac{1}{2}}h^{-\frac{1}{2}}\|\tilde{e}_u\| \\ &\le c\|\tilde{e}_u\|^2 + ch^{2k}. \end{aligned} \quad (54)$$

Thus by Gronwall's inequality, the conclusion in Theorem B.2 follows. □

## 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 = \underbrace{[1, \underbrace{0, \dots, 0}_{\in \mathbb{R}^k}, \underbrace{1, \underbrace{0, \dots, 0}_{\in \mathbb{R}^k}, \underbrace{\dots, \underbrace{1, \underbrace{0, \dots, 0}_{\in \mathbb{R}^k}}, \dots]}_{\in \mathbb{R}^{(k+1)N}}]^T \in \mathbb{R}^{(k+1)N}.$$

is an 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, \mathbf{0}, \dots, \mathbf{0}, u_1, \mathbf{0}, \dots, \mathbf{0}, \dots, u_1, \mathbf{0}, \dots, \mathbf{0}] \in \mathbb{M}_{(k+1)N, (k+1)N}(\mathbb{R}),$$

with $\mathbf{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 = \underbrace{[0, 1, 0, 1, \dots, 0, 1]}_{\in \mathbb{R}^2}^T \in \mathbb{R}^N.$$

$$k=1, u_2 = \underbrace{[\underbrace{0, 1}_{\in \mathbb{R}^2}, \underbrace{0, 1}_{\in \mathbb{R}^2}, \dots, \underbrace{0, 1}_{\in \mathbb{R}^2}]}_{\in \mathbb{R}^{2N}}^T \in \mathbb{R}^{2N}.$$
$$k = 2, u_2 = \left[ \underbrace{-\frac{1}{6}, 0, 1, 0, 0, -1}_{\in \mathbb{R}^6}, \dots, \underbrace{-\frac{1}{6}, 0, 1, 0, 0, -1}_{\in \mathbb{R}^6} \right]^T \in \mathbb{R}^{3N}.$$

$$k = 3, u_2 = \left[ \underbrace{0, -\frac{3}{20}, 0, 1, 0, -\frac{3}{20}, 0, 1, \dots, 0, -\frac{3}{20}, 0, 1}_{\in \mathbb{R}^4} \right]^T \in \mathbb{R}^{4N}.$$

$$k = 4, u_2 = \left[ \underbrace{-\frac{3}{280}, 0, \frac{3}{14}, 0, -1, 0, 0, -\frac{3}{14}, 0, 1, \dots, \underbrace{-\frac{3}{280}, 0, \frac{3}{14}, 0, -1, 0, 0, -\frac{3}{14}, 0, 1}_{\in \mathbb{R}^{10}} \right]^T \in \mathbb{R}^{5N}.$$

$$k = 5, u_2 = \left[ \underbrace{0, \frac{5}{336}, 0, -\frac{5}{18}, 0, 1, 0, \underbrace{\frac{5}{336}, 0, -\frac{5}{18}, 0, 1, \dots, 0, \frac{5}{336}, 0, -\frac{5}{18}, 0, 1}_{\in \mathbb{R}^6} \right]^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$.

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

In this appendix, we extend the method presented in [5] 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 (24) satisfied by $f(t, x, v_1, v_2)$, $E_1(t, x)$, $E_2(t, x)$, $B(t, x)$, with $t \ge 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} i k \hat{f}_{k,j_1,j_2} + (\widehat{\mathcal{F}Df})_{k,j_1,j_2} &= 0, && \text{with } \mathcal{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.
\qquad (55)$$

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 \ne 0$. Let us denote $U = (\hat{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), \qquad (56)$$
with

$$
\begin{align*}
& U = (\hat{f}_{k,\star,1}, \hat{f}_{k,\star,2}, \dots, \hat{f}_{k,\star,N_{v_2}}, \hat{B}_k, \hat{E}_{2,k}, \hat{E}_{1,k})^T, \\
& \hat{f}_{k,\star,j_2} \in \mathbb{C}^{N_{v_1}}, \forall k = 1, \dots, N_x, \forall j_2 = 1, \dots, N_{v_2}, \\
& N(U) = (-\widehat{(\mathcal{FD}f)}_{k,\star,1}, -\widehat{(\mathcal{FD}f)}_{k,\star,2}, \dots, -\widehat{(\mathcal{FD}f)}_{k,\star,N_{v_2}}, 0, 0, 0)^T, \\
& L = \begin{pmatrix}
-ik\operatorname{diag}(\vec{v_1}) & \mathbf{0}_{N_{v_1},N_{v_1}} & \mathbf{0}_{N_{v_1},N_{v_1}} & \dots & \mathbf{0}_{N_{v_1},N_{v_1}} & \mathbf{0}_{N_{v_1},1} & \dots & \mathbf{0}_{N_{v_1},1} \\
\mathbf{0}_{N_{v_1},N_{v_1}} & -ik\operatorname{diag}(\vec{v_1}) & \mathbf{0}_{N_{v_1},N_{v_1}} & \dots & \mathbf{0}_{N_{v_1},N_{v_1}} & \vdots & \vdots & \vdots \\
\vdots & \ddots & \vdots & \ddots & \mathbf{0}_{N_{v_1},N_{v_1}} & \vdots & \vdots & \vdots \\
\mathbf{0}_{N_{v_1},N_{v_1}} & \dots & \mathbf{0}_{N_{v_1},N_{v_1}} & -ik\operatorname{diag}(\vec{v_1}) & \mathbf{0}_{N_{v_1},1} & \vdots & \mathbf{0}_{N_{v_1},1} & \vdots \\
\mathbf{0}_{1,N_{v_1}} & \dots & \mathbf{0}_{1,N_{v_1}} & \mathbf{0}_{1,N_{v_1}} & \mathbf{0}_{1,N_{v_1}} & 0 & -ik & 0 \\
-\Delta v_1 \Delta v_2 (\vec{v_2})_1 & \mathbf{1} & -\Delta v_1 \Delta v_2 (\vec{v_2})_2 & \mathbf{1} & \dots & -\Delta v_1 \Delta v_2 (\vec{v_2})_{N_{v_2}} & \mathbf{1} & -ik & 0 \\
-\Delta v_1 \Delta v_2 \vec{v_1} & \mathbf{1} & -\Delta v_1 \Delta v_2 \vec{v_1} & \mathbf{1} & \dots & -\Delta v_1 \Delta v_2 \vec{v_1} & 0 & 0 & 0
\end{pmatrix} \tag{57}
\end{align*}
$$

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, $\operatorname{diag}(\vec{v}_1)$ denotes the diagonal matrix whose diagonal is composed of $\vec{v}_1$, $\mathbf{1} \in \mathbb{R}^{N_{v_1}}$ denotes the vector with components 1 and $\mathbf{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 (56). 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 (57) is given by

$$
\begin{pmatrix}
e^{-ik\Delta t \operatorname{diag}(\vec{v}_1)} & \mathbf{0}_{N_{v_1},N_{v_1}} & \mathbf{0}_{N_{v_1},N_{v_1}} & \dots & \mathbf{0}_{N_{v_1},N_{v_1}} & \mathbf{0}_{N_{v_1},1} & \dots & \mathbf{0}_{N_{v_1},1} \\
\mathbf{0}_{N_{v_1},N_{v_1}} & e^{-ik\Delta t \operatorname{diag}(\vec{v}_1)} & \mathbf{0}_{N_{v_1},N_{v_1}} & \dots & \mathbf{0}_{N_{v_1},N_{v_1}} & \vdots & \vdots & \vdots \\
\vdots & \ddots & \vdots & \ddots & \mathbf{0}_{N_{v_1},N_{v_1}} & \vdots & \vdots & \vdots \\
\mathbf{0}_{N_{v_1},N_{v_1}} & \dots & \mathbf{0}_{N_{v_1},N_{v_1}} & e^{-ik\Delta t \operatorname{diag}(\vec{v}_1)} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} & \mathbf{0}_{N_{v_1},1} \\
e \mathcal{B}_{k,1} & e \mathcal{B}_{k,2} & \dots & e \mathcal{B}_{k,N_{v_2}} & \cos(k\Delta t) & -i \sin(k\Delta t) & 0 & 0 \\
e \mathcal{E}_{2,k,1} & e \mathcal{E}_{2,k,2} & \dots & e \mathcal{E}_{2,k,N_{v_2}} & -i \sin(k\Delta t) & \cos(k\Delta t) & 0 & 0 \\
e \mathcal{E}_{1,k} & e \mathcal{E}_{1,k} & \dots & e \mathcal{E}_{1,k} & 0 & 0 & 1 & 1
\end{pmatrix}
$$

where $e \mathcal{B}_{k,j_2}, e \mathcal{E}_{2,k,j_2}, e \mathcal{E}_{1,k} \in \mathbb{M}_{1,N_{v_1}}(\mathbb{C})$ for $j_2 = 1, \dots, N_{v_2}$ are given by

$$e \mathcal{B}_{k,j_2} = \frac{i \Delta v_1 \Delta v_2}{k} v_{2,j_2} \vec{\alpha}, \quad e \mathcal{E}_{2,k,j_2} = \frac{i \Delta v_1 \Delta v_2}{k} v_{2,j_2} \vec{\beta}, \quad e \mathcal{E}_{1,k} = \frac{i \Delta v_1 \Delta v_2}{k} (\mathbf{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

$$
\begin{align*}
\alpha_{j_1} &= -\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 t v_{1,j_1}}}{(1 - v_{1,j_1}^2)}, \\
\beta_{j_1} &= -\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 t v_{1,j_1}}}{(1 - v_{1,j_1}^2)}.
\end{align*}
$$
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\bar{v}_1}$ ($N_{v_2}$ times). Second, the 3 × 3 right bottom block corresponds to the homogeneous Maxwell equations and its exponential can be computed as

$$ \exp\begin{pmatrix} 0 & -ik\Delta t & 0 \\ -ik\Delta t & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} \cos(k\Delta t) & i\sin(k\Delta t) & 0 \\ i\sin(k\Delta t) & \cos(k\Delta t) & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$

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

**Computation of $e^{i\mathcal{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})$

$$ \begin{aligned} \hat{E}_{1,k}(t^{n+1}) &= \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)} dt v_{1,j_1} f_{k,j_1,j_2}(t^n) \Delta v_1 \Delta v_2 \\ &= \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), \end{aligned} $$

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

**Computation of $e^{i\mathcal{B}_k}$, $e^{i\mathcal{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

$$ \frac{d}{dt} \hat{E}_{2,k}(t) = -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 $$

$$ \frac{d}{dt} \hat{B}_k(t) = -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 = \begin{pmatrix} 0 & -ik \\ -ik & 0 \end{pmatrix} \quad \text{and} \quad R = \begin{pmatrix} R_1(t) \\ 0 \end{pmatrix}, \quad \text{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. \quad (58) $$

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

$$ e^{M\Delta t} = \exp\begin{pmatrix} 0 & -ik\Delta t \\ -ik\Delta t & 0 \end{pmatrix} = \begin{pmatrix} \cos(kt) & -i\sin(kt) \\ -i\sin(kt) & \cos(kt) \end{pmatrix}. $$
Second, one has to compute the integral term in (58)

$$
\begin{align}
\int_{t^n}^{t^{n+1}} e^{-M(t-t^{n+1})} R(t) dt &= \left( \begin{array}{c} \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) \nonumber \\
&= - \left( \begin{array}{c} \sum_{j_1, j_2} \mathcal{I}_1 v_{2,j_1} \hat{f}_{k,j_1,j_2}^n \Delta v_1 \Delta v_2 \\ \sum_{j_1, j_2} \mathcal{I}_2 v_{2,j_1} \hat{f}_{k,j_1,j_2}^n \Delta v_1 \Delta v_2 \end{array} \right) \tag{59}
\end{align}
$$

where $\mathcal{I}_1$ and $\mathcal{I}_2$ are given by

$$
\mathcal{I}_1 = \int_{t^n}^{t^{n+1}} \left[ \cos(k(t - t^{n+1})) e^{-ikv_{1,j_1}(t-t^n)} \right] 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)},
$$

$$
\mathcal{I}_2 = \int_{t^n}^{t^{n+1}} \left[ i \sin(k(t - t^{n+1})) e^{-ikv_{1,j_1}(t-t^n)} \right] 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 (58) leads to the following expression for $\hat{E}_{2,k}(t)$

$$
\hat{E}_{2,k}(t^{n+1}) = \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),
$$

$$
\hat{B}_k(t^{n+1}) = -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

$$
\beta_{j_1} = - \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)},
$$

$$
\alpha_{j_1} = - \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\mathcal{B}_k, e\mathcal{E}_{2,k}$ corresponding to $\hat{B}_k$ and $\hat{E}_{2,k}$

$$
e\mathcal{B}_{k,j_2} = \frac{i\Delta v_1 \Delta v_2}{k} v_{2,j_2} \vec{\alpha}, \quad e\mathcal{E}_{2,k,j_2} = \frac{i\Delta v_1 \Delta v_2}{k} v_{2,j_2} \vec{\beta}.
$$

□
References

[1] B. Ayuso de Dios, J. A. Carrillo de la Plata, and C.-W. Shu. Discontinuous Galerkin methods for the one-dimensional Vlasov-Poisson system. 2009.

[2] 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.

[3] 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.

[4] 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.

[5] B. Boutin, A. Crestetto, N. Crouseilles, and J. Massot. Modified lawson methods for Vlasov equations. 2022.

[6] 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.

[7] 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.

[8] 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.

[9] 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.

[10] Y. Cheng, A. J. Christlieb, and X. Zhong. Energy-conserving discontinuous Galerkin methods for the Vlasov-Ampere system. *Journal of Computational Physics*, 256:630-655, 2014.

[11] 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.

[12] Y. Cheng, I. M. Gamba, A. Majorana, and C.-W. Shu. Discontinuous Galerkin solver for Boltzmann-Poisson transients. *Journal of Computational Electronics*, 7:119-123, 2008.

[13] P. G. Ciarlet. *The finite element method for elliptic problems*. SIAM, 2002.

[14] B. Cockburn, F. Li, and C.-W. Shu. Locally divergence-free discontinuous Galerkin methods for the Maxwell equations. *Journal of Computational Physics*, 194(2):588-610, 2004.

[15] B. Cockburn and C.-W. Shu. Runge-Kutta discontinuous Galerkin methods for convection-dominated problems. *Journal of scientific computing*, 16:173-261, 2001.

[16] 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.
[17] 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.

[18] 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.

[19] B. A. de Dios and S. Hajian. High order and energy preserving discontinuous Galerkin methods for the Vlasov-Poisson system. *arXiv preprint arXiv:1209.4025*, 2012.

[20] 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.

[21] F. Filbet and E. Sonnendrücker. Comparison of Eulerian Vlasov solvers. *Computer Physics Communications*, 150(3):247-266, 2003.

[22] M. Hochbruck and A. Ostermann. Exponential integrators. *Acta Numerica*, 19:209-286, 2010.

[23] X. Hong and J.-M. Qiu. A generalized Eulerian-Lagrangian discontinuous Galerkin method for transport problems. *Journal of Computational Physics*, 464:111160, 2022.

[24] 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.

[25] 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.

[26] 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.

[27] 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.

[28] 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.

[29] 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.

[30] 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.

[31] J.-M. Qiu and A. Christlieb. A conservative high order semi-Lagrangian WENO method for the Vlasov equation. *Journal of Computational Physics*, 229(4):1130-1149, 2010.

[32] G. R. Richter. An optimal-order error estimate for the discontinuous Galerkin method. *Mathematics of Computation*, 50(181):75-88, 1988.
[33] J. A. Rossmanith and D. C. Seal. A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations. *Journal of Computational Physics*, 230(16):6203-6232, 2011.

[34] G. J. Tee. Eigenvectors of block circulant and alternating circulant matrices. *New Zealand Journal of Mathematics*, 36(8):195-211, 2007.

[35] J. P. Verboncoeur. Particle simulation of plasmas: review and advances. *Plasma Physics and Controlled Fusion*, 47(5A):A231, 2005.

[36] H. Wang, R. E. Ewing, G. Qin, S. L. Lyons, M. Al-Lawatia, and S. Man. A family of eulerian-lagrangian localized adjoint methods for multi-dimensional advection-reaction equations. *Journal of Computational Physics*, 152(1):120-163, 1999.

[37] E. S. Weibel. Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution. *Physical Review Letters*, 2(3):83, 1959.

[38] Q. Zhang and C.-W. Shu. Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws. *SIAM Journal on Numerical Analysis*, 42(2):641-666, 2004.

[39] Q. Zhang and C.-W. Shu. Stability analysis and a priori error estimates of the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws. *SIAM Journal on Numerical Analysis*, 48(3):1038-1063, 2010.