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

HAL
open science

# Modified lawson methods for vlasov equations *

Benjamin Boutin, Anais Crestetto, Nicolas Crouseilles, Josselin Massot

► To cite this version:

Benjamin Boutin, Anais Crestetto, Nicolas Crouseilles, Josselin Massot. Modified lawson methods for vlasov equations *. 2022. hal-03911409

HAL Id: hal-03911409

https://hal.science/hal-03911409

Preprint submitted on 22 Dec 2022

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.
MODIFIED LAWSON METHODS FOR VLASOV EQUATIONS*

BENJAMIN BOUTIN†, ANAÏS CRESTETTO‡, NICOLAS CROUSEILLES§, AND JOSSELIN MASSOT¶

**Abstract.** In this work, Lawson type numerical methods are studied to solve Vlasov type equations on a phase space grid. These time integrators are known to satisfy enhanced stability properties in this context since they do not suffer from the stability condition induced from the linear part. We introduce here a class of modified Lawson integrators in which the linear part is approximated in such a way that some geometric properties of the underlying model are preserved, which has important consequences for the analysis of the scheme. Several Vlasov-Maxwell examples are presented to illustrate the good behavior of the approach.

**Key words.** Lawson methods, Vlasov equations, high order.

**MSC codes.** 35L45, 65L06, 65L07, 65M12

**1. Introduction.** In this work, we study high order Eulerian numerical methods to solve Vlasov type equations describing a system of charged particles under the influence of self-consistent electromagnetic fields. In addition to their nonlinear character, the solutions of the Vlasov equations present fine structures and strong gradients which are difficult to capture for the typical time scales used in plasma physics.

Several methods have been proposed in the past to solve efficiently kinetic problems. First, the so-called Particle-In-Cell methods have been introduced and still are extensively used [7, 48, 34, 33]. They are based on a sampling of the distribution function using macro particles which are advanced in time through the characteristics equations, whereas the electromagnetic fields are computed on a spatial grid thanks to grid-projection techniques of the macro-particles. Even if these methods are extremely efficient in high dimensions, they are known to suffer from some numerical noise which prevents from an accurate description of the unknown and which slowly decreases when the number of macro-particles increases. Second, the so-called Eulerian methods which use a grid of the phase space have been considered [10, 32, 25, 44, 2, 42, 43, 23, 49, 13, 9, 3, 15, 6]. These approaches enable the use of high order methods like spectral or finite differences/volumes methods, and as such, capture fine physical phenomena like Landau damping or filamentation. However, the stability condition relating the phase space mesh and the time step might be stringent, making these methods costly both in terms of CPU and memory to reach typical physical time scales. Another approach called the semi-Lagrangian method has been considered to take the best of the two latter approaches ([45, 25, 18, 17]). This method is Eulerian since it requires a grid of the phase space but it is Lagrangian since it follows the characteristics equation to update the unknown from one iteration to the next as such, the method does not suffer from the stability condition. How-

<footnote>*Submitted to the editors January 2023.</footnote>

<footnote>†Univ Rennes, CNRS, IRMAR - UMR 6625, F-35000 Rennes, France (benjamin.boutin@universennes1.fr).</footnote>

<footnote>‡Nantes Université, CNRS, Laboratoire de Mathématiques Jean Leray, LMJL, UMR 6629, F-44000 Nantes, France (anais.crestetto@univ-nantes.fr).</footnote>

<footnote>§Univ Rennes, IRMAR UMR 6625 & centre Inria de l'Université de Rennes (MINGuS) & ENS Rennes, France (nicolas.crouseilles@inria.fr).</footnote>

<footnote>¶CMAP, CNRS, Ecole polytechnique, Institut Polytechnique de Paris, France (josselin.massot@polytechnique.edu)</footnote>
ever, to avoid high dimensional costly phase space reconstructions or interpolations, they are very often combined with a time splitting. Splitting methods turn out to be very powerful for Vlasov problems since the subsystems can often be solved exactly in time. However, it is not always the case and high order in time requires a very large number of stages, which of course may lead to expensive methods, in particular when the splitting involves three or more subsystems (see [12]).

Recently (see [16, 15, 12, 39]), exponential time integrators have been investigated to numerically solve Vlasov equations as an alternative to above methods. The main observation was based on the fact that, in several Vlasov equations, the linear part induces the most stringent CFL condition. But since, by essence, the linear part can be potentially solved easily, a variation of constant formula can be written to pave the way of a new class of methods where the linear part is solved exactly and the nonlinear part is solved explicitly in time (to avoid expensive iterative solvers). Thanks to the huge bibliography on Runge-Kutta (RK) methods, efficient high order methods can then be obtained, whose number of stages (and thus the complexity) is only linear with respect to the order of the time integrator. There exists two main families of exponential integrators (see [29]) and as investigated in [15], among them, Lawson methods [36] are preferred for Vlasov problems.

In this work, we investigate and analyse Lawson methods for different Vlasov type equations. More specifically, for some complex problems, the exponential of the linear part turns out to be difficult (even impossible) to compute, even using dedicated softwares and some approximations are required. But due to the specific properties of the linear part in Vlasov equations, some standard exponential approximations lead to unstable methods. Indeed, when one considers the Vlasov equations, the linear part has pure imaginary eigenvalues, which means that the eigenvalues of the exponential belong to the unit circle, for any time steps and any space meshes. This important geometric property has to be preserved by the approximation for stability reasons. We checked that the so-called Padé approximant satisfies this property and the Padé strategy turns out to be well adapted in our context. Then, we prove that the order of accuracy of the so modified Lawson method (in which the exact exponential is replaced by an approximation) becomes $O(\Delta t^{\min(m,r)})$ with $\Delta t$ the time step, $m$ the order of the Lawson method and $r$ the order of approximation of the exponential approximation. As these methods rely on the very popular RK methods, they benefit from the huge literature enabling to reach high order efficiently and to be combined with adaptive time stepping strategy by using embedded RK methods.

In addition, for more simple problems like Vlasov-Ampère or Vlasov-Maxwell, for which the exponential of the linear part can be computed exactly, we observe that Lawson schemes enable to preserve the charge conservation exactly, for any underlying RK method. Indeed, inserting the calculation of the current in the linear part enables to solve exactly both the space transport and the Ampère equation, which is the key point to ensure that the underlying Poisson equation is satisfied exactly (see [44, 33, 19]). Hence, in addition to its efficiency, the versatility of exponential integrators turns out to be an interesting aspect in our context.

The rest of the paper is organized as follows. First, several aspects of Lawson methods are recalled and modified Lawson methods (where the exponential of the matrix is approximated) is introduced. Second, some models amenable to Lawson integrators are presented and in a third part, the Padé approximants are recalled and some important properties in our context are proved. Then, we analyse the convergence of the modified Lawson methods for a class of ODEs. Finally, several numerical illustrations are discussed.
**2. Lawson methods.** In this section, we recall the Lawson methods and introduce the modified Lawson methods. Lawson methods are a class of time integration schemes that are applied to differential equations (which, in our context, come from the phase space semi-discretization of Vlasov equations) of the form:

$$ (2.1) \qquad \dot{u}(t) = Lu(t) + g(t, u(t)), \quad u(0) = u_0 \in \mathbb{R}^d, $$

where $L$ is a $d \times d$ matrix and $g$ is a, in general nonlinear, function of the unknown $u$ and the time $t \ge 0$ whose values belong to $\mathbb{R}^d$ (i.e. for $d \in \mathbb{N}$, we have $u : t \in \mathbb{R}^+ \mapsto u(t) \in \mathbb{R}^d$ and $g : (t, u) \in \mathbb{R}^+ \times \mathbb{R}^d \mapsto \mathbb{R}^d$). Lawson methods are especially efficient when applied to problems where $L$ implies a stringent stability condition if it is treated explicitly, which is typically the case of the Vlasov equations (see [15, 12]).

**2.1. Standard Lawson methods.** Solving (2.1) with the Lawson method [35] relies on the change of variable $v(t) := \exp(-tL)u(t)$ in (2.1) to get

$$ (2.2) \qquad \dot{v}(t) = \exp(-tL)g(t, \exp(tL)v(t)). $$

Now an explicit Runge-Kutta method is applied to the transformed equation (2.2). We introduce the time discretization $t^n = n\Delta t$ with $\Delta t > 0$ the time step, $n \in \mathbb{N}$ and $u^n$ (resp. $v^n$) denotes the numerical approximation of $u(t^n)$ (resp. $v(t^n)$). For the sake of simplicity, we first present the method for the explicit Euler scheme. Applying the forward Euler method to (2.2) leads to

$$ v(t^n + \Delta t) \approx v^{n+1} = v^n + \Delta t \exp(-t^n L) g(t^n, \exp(t^n L) v^n). $$

Reversing the change of variable yields the following scheme for $u^n$

$$ u(t^n + \Delta t) \approx u^{n+1} = \exp(\Delta t L) u^n + \Delta t \exp(\Delta t L) g(t^n, u^n). $$

This is the Lawson-Euler method, also a method of order one. More generally, the Lawson method induced by an explicit Runge-Kutta method RK(s, m) with s stages of order m can be written as (2.3)

$$ u_n^{(r)} = \exp(c_r \Delta t L) u^n + \Delta t \sum_{j=1}^{r-1} a_{r,j} \exp((c_r - c_j) \Delta t L) g(t^n + c_j \Delta t, u_n^{(j)}), \quad r = 1, \dots, s, $$

$$ u^{n+1} = \exp(\Delta t L) u^n + \Delta t \sum_{j=1}^{s} b_j \exp((1 - c_j) \Delta t L) g(t^n + c_j \Delta t, u_n^{(j)}), $$

with the underlying explicit Runge-Kutta method defined by its Butcher tableau. The method defined by (2.3) will be denoted by LRK(s, m), which is the Lawson method induced by the RK(s, m) Runge-Kutta method.

**2.2. Modified Lawson methods.** For many problems, $\exp(tL)$ may be difficult, costly or even impossible to compute formally for all the involved parameters and approximations may be required (see [8, 22, 24, 40]). In the sequel we will denote by $\exp(p(tL))$ an approximation of $\exp(tL)$ and introduce $\tilde{u}^n$ an approximation of $u(t^n)$ for $n \in \mathbb{N}$ generated by the modified Lawson method. The modified Lawson M-LRK(s, m) method (of order m with s stages) is simply obtained from (2.3) where
the matrix exponentials $\exp(tL)$ is approximated by its approximation $\mathfrak{erp}(tL)$ to get (2.4)

$$ \tilde{u}_n^{(r)} = \mathfrak{erp}(c_r \Delta t L) \tilde{u}^n + \Delta t \sum_{j=1}^{r-1} a_{r,j} \mathfrak{erp}((c_r - c_j) \Delta t L) g(t^n + c_j \Delta t, \tilde{u}_n^{(j)}), \quad r = 1, \dots, s, $$

$$ \tilde{u}^{n+1} = \mathfrak{erp}(\Delta t L) \tilde{u}^n + \Delta t \sum_{j=1}^{s} b_j \mathfrak{erp}((1 - c_j) \Delta t L) g(t^n + c_j \Delta t, \tilde{u}_n^{(j)}), $$

where the coefficients are given by the Butcher tableau.

**3. Some Vlasov models amenable to Lawson methods.** In this section, we present some Vlasov type models that enter in the framework presented in the previous section. We focus here on Eulerian methods which consider a grid of the phase space to solve the Vlasov equations. For the different Vlasov models we consider, a reformulation into the form (2.1) is proposed so that Lawson or modified Lawson methods can be employed. We discuss the eigenvalues of the linear part $L$. Indeed, in our case $L$ has a special structure since it has a pure imaginary spectrum which has important consequences for the stability of the Lawson schemes studied here.

**3.1. Vlasov-Ampère.** The Vlasov-Ampère equation for the distribution function $f(t, x, v) \ge 0$ and the electric field $E(t, x) \in \mathbb{R}$ with $x \in [0, b]$ and $v \in \mathbb{R}$ is

$$ \frac{\partial f}{\partial t} + v \partial_x f + E \partial_v f = 0, \quad \frac{\partial E}{\partial t} = - \int_{\mathbb{R}} v f dv + \frac{1}{b} \int_0^b \int_{\mathbb{R}} v f dv dx, $$

with the initial conditions ($f_0, E_0$) such that the Poisson equation is satisfied initially (i.e. $\partial_x E_0 = \int_{\mathbb{R}} f_0 dv - 1$) and periodic boundary conditions are imposed in space.

Let us detail the numerical strategy in this case. We first semi-discretize the phase space domain $(x, v)$: we have chosen spectral method in space and high order finite differences in velocity. The velocity mesh is $v_\ell = -v_{\max} + (\ell - 1)\Delta v$, $\ell = 1, \dots, N_v$ where $v_{\max} > 0$ is the velocity domain truncation and $\Delta v = (2v_{\max})/(N_v - 1)$ is the velocity mesh size. Then, we denote $\tilde{f}_k(t) = (\hat{f}_{k,1}, \hat{f}_{k,2}, \dots, \hat{f}_{k,N_v})(t)$ the semi-discrete unknown for $f$ and $\hat{E}_k(t)$ the semi-discrete unknown for $E$, $k = 0, \dots, N_x - 1$ ($N_x$ being the number of points of the space mesh): (3.1)

$$ \hat{f}_{k,\ell}(t) = \frac{1}{b} \int_0^b \exp(-i2\pi kx/b) f(t, x, v_\ell) dx, \quad \hat{E}_k(t) = \frac{1}{b} \int_0^b \exp(-i2\pi kx/b) E(t, x) dx. $$

We will introduce the following notations: $\vec{v} = (v_1, v_2, \dots, v_{N_v}) \in \mathbb{R}^{N_v}$ for the mesh in velocity and $\mathbf{0} = (0, 0, \dots, 0) \in \mathbb{R}^{N_v}$. We then get

$$ (3.2) \quad \frac{d}{dt} \begin{pmatrix} \vec{f}_k(t) \\ \hat{E}_k(t) \end{pmatrix} + \begin{pmatrix} ik \frac{2\pi}{b} \operatorname{diag}(\vec{v}) & \mathbf{0}^T \\ \Delta v \vec{v} & 0 \end{pmatrix} \begin{pmatrix} \vec{f}_k(t) \\ \hat{E}_k(t) \end{pmatrix} + \left( \left[ \widehat{(ED\vec{f})} \right]_k \right) = 0, $$

with $\operatorname{diag}(\vec{v})$ is the diagonal matrix of size $N_v$ with $v_1, v_2, \dots, v_{N_v}$ on its diagonal and $(EDf) \in \mathbb{R}^{N_v}$ is a finite difference approximation of $E\partial_v f$ on the velocity mesh. The (Fourier coefficient of the) current in the Ampère equation is approximated by a standard quadrature $\sum_{\ell=1}^{N_v} v_\ell \hat{f}_{k,\ell}(t)\Delta v$. We can then reformulate (3.2) under the following form with $U_k(t) = (\tilde{f}_k(t), \hat{E}_k(t)) \in \mathbb{R}^{N_v+1}$

$$ \frac{d}{dt} U_k(t) = LU_k(t) + N(U_k(t)), $$
139 which is amenable to Lawson schemes. Moreover, $L$ is diagonalizable $L = PDP^{-1}$
140 with a pure imaginary spectrum $\text{Sp}(L) = \{i(2\pi k/b)\vec{v}\}$ and $\text{cond}(P)$ is uniformly
141 bounded with respect to $N_v$ and $N_x$.

142 Let us remark that this approach is different from the one proposed in [15, 16]
143 since here the Ampère equation is involved in the linear part whereas the Poisson
144 equation was updated at each Runge-Kutta stage in [15, 16]. We believe that the
145 current approach is more amenable to the Vlasov-Maxwell equation; moreover, as we
146 shall see in Section 6, some properties can be proved using this approach.

**3.2. Linearized hybrid fluid-kinetic model.** The second example considered here is a hybrid model introduced in [46, 47, 30]. In this model, the case of a wave propagation parallel to a uniform magnetic field $\mathbf{B}_0 = (0,0, B_0)^T$, $B_0 > 0$ is considered so that the problem becomes one dimensional in space (called $z$) but the three dimensions in velocity are kept: $\mathbf{v} = (v_x, v_y, v_z) \in \mathbb{R}^3$. Thus, under some assumptions (see [30, 12]), a hot/cold decomposition is performed, and after a cold plasma approximation for the cold electron population, the following decomposition of the distribution function is considered:

$$f(t, z, \mathbf{v}) = \rho_c(z)\delta(\mathbf{v} - \mathbf{u}_c(t, z)) + f_h(t, z, \mathbf{v}),$$

with $\mathbf{j}_c(t, z) = \rho_c(z)\mathbf{u}_c(t, z) = (j_{c,x}, j_{c,y}, 0)(t, z) \in \mathbb{R}^3$ denotes the current, $\mathbf{u}_c(t, z) \in \mathbb{R}^3$ the mean velocity and $\rho_c(z) \in \mathbb{R}$ the density of the cold particles population. Hence, the linearized hybrid fluid-kinetic model satisfied by

$$U(t, z, \mathbf{v}) = (f_h(t, z, \mathbf{v}), E_x(t, z), E_y(t, z), B_x(t, z), B_y(t, z), j_{c,x}(t, z), j_{c,y}(t, z)),$$

147 can be derived (see [30, 12] for more details). The model is normalized as in [30, 12]
148 and reads as $(\Omega_{pe})$ denotes the ratio between the plasma and cyclotronic frequencies)

$$149 \quad (3.3) \quad \frac{\partial j_{c,x}}{\partial t} = \Omega_{pe}^2 E_x - j_{c,y} B_0, \quad \frac{\partial j_{c,y}}{\partial t} = \Omega_{pe}^2 E_y + j_{c,x} B_0,$$

$$150 \quad (3.4) \quad \frac{\partial B_x}{\partial t} = \partial_z E_y, \quad \frac{\partial B_y}{\partial t} = -\partial_z E_x,$$

$$151 \quad (3.5) \quad \frac{\partial E_x}{\partial t} = -\partial_z B_y - j_{c,x} + \int_{\mathbb{R}^3} v_x f_h \, d\mathbf{v}, \quad \frac{\partial E_y}{\partial t} = \partial_z B_x - j_{c,y} + \int_{\mathbb{R}^3} v_y f_h \, d\mathbf{v},$$

$$152 \quad (3.6) \quad \frac{\partial f_h}{\partial t} + v_z \partial_z f_h - \mathcal{F} \cdot \nabla_{\mathbf{v}} f_h = 0,$$

$$154 \quad \text{with } \mathcal{F}(t, z, v_x, v_y, v_z) = (E_x + v_y B_0 - v_z B_y, E_y - v_x B_0 + v_z B_x, v_x B_y - v_y B_x).$$

155 As in the previous case, we use a Fourier discretization in the space $z$ direction
156 so that the unknown is $U_k(t) = (\hat{j}_{c,x}, \hat{j}_{c,y}, \hat{B}_x, \hat{B}_y, \hat{E}_x, \hat{E}_y, \hat{f}_h)_k$ and satisfies $\partial_t U_k =$
157 $LU_k + N(t, U_k)$ with
(3.7)

$$158 \quad L = \begin{pmatrix} 0 & -B_0 & 0 & 0 & \Omega_{pe}^2 & 0 & 0 \\ B_0 & 0 & 0 & 0 & 0 & \Omega_{pe}^2 & 0 \\ 0 & 0 & 0 & 0 & 0 & ik & 0 \\ 0 & 0 & 0 & 0 & -ik & 0 & 0 \\ -1 & 0 & 0 & -ik & 0 & 0 & 0 \\ 0 & -1 & ik & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -v_z ik \end{pmatrix}, \quad N: (t,U) \mapsto \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ \int_{\mathbb{R}^3} v_x (\hat{f}_h)_k \, d\mathbf{v} \\ \int_{\mathbb{R}^3} v_y (\hat{f}_h)_k \, d\mathbf{v} \\ [\mathcal{F} \cdot \nabla_{\mathbf{v}} f_h]_k \end{pmatrix}.$$
For the velocity space, we use finite differences methods (WENO type). Moreover, $L$ is diagonalizable $L = PDP^{-1}$ with a pure imaginary spectrum and $\text{cond}(P)$ is uniformly bounded with respect to $N_x, N_v$.

**4. Padé approximant.** In this section we recall the Padé approximant to approximate exponential of matrices. Obviously, many other approximations are possible (see [41]) but the Padé approximant enjoys some properties that have important consequences for the numerical approximation of the Vlasov models presented in the previous section in terms of stability of the resulting numerical schemes.

The Padé method turns out to be a very popular technique to approximate exponential of matrices (see for instance [1, 38, 4]). For a given matrix $L$, $t \in \mathbb{R}$ and $(p, q) \in \mathbb{N}^2$, the Padé approximant is a rational function defined by

$$ (4.1) \qquad P_{p,q}(tL) = h_{p,q}(tL) (k_{p,q}(tL))^{-1}, $$

the numerator and denominator being polynomial functions of degree $p$ and $q$

$$ h_{p,q}(tL) = \sum_{\ell=0}^{p} \frac{\frac{p!}{(p-\ell)!}}{\frac{(p+q)!}{(p+q-\ell)!}} \frac{(tL)^{\ell}}{\ell!}, \qquad k_{p,q}(tL) = \sum_{\ell=0}^{q} (-1)^{\ell} \frac{\frac{q!}{(q-\ell)!}}{\frac{(p+q)!}{(p+q-\ell)!}} \frac{(tL)^{\ell}}{\ell!}. $$

From [1, 41], the error between the exponential of the matrix and the Padé approximant satisfies (for a given matrix norm $\| \cdot \|$)

$$ (4.2) \qquad \| \exp(tL) - P_{p,q}(tL) \| \le C t^{p+q+1}. $$

The choice of the Padé approximant is motivated in our case by the following proposition.

**PROPOSITION 4.1.** Let $L$ be a diagonalizable matrix of size $d \in \mathbb{N}^*$ (there exist an invertible matrix $P$ and a diagonal matrix $D$ such that $L = PDP^{-1}$) such that $\text{Sp}(L) \subset i\mathbb{R}$ and let assume there exists a constant $C > 0$ such that $\text{cond}(P) \le C$. Then, for $p \in \mathbb{N}$, the eigenvalues of the Padé approximant $P_{p,p}(tL)$ of $\exp(tL)$ belongs to the unit circle $C(0,1)$, i.e., $\forall \lambda_j \in \text{Sp}(P_{p,p}(tL))$, one has $|\lambda_j| = 1$. Moreover, for any $t \in \mathbb{R}$ and for any $n \in \mathbb{N}$

$$ (4.3) \qquad \|(P_{p,p}(tL))^n\| \le C. $$

*Proof.* Because $L$ is diagonalizable, there exists an invertible matrix $P$ such that $L = PDP^{-1}$ with $D$ a diagonal matrix. Since $\text{Sp}(L) \subset i\mathbb{R}$, the diagonal terms of $D$ can be written as $D_{j,j} = i\alpha_j$, $\alpha_j \in \mathbb{R}$, $j = 1, \dots, d$. Computing powers of $L$ thus reduces to look at the power of $D_{j,j}$.

When $p = q$ the Padé approximant $P_{p,p}(tL)$ can be written as

$$ P_{p,p}(tL) = \left( \sum_{\ell=0}^{p} a_{\ell,p} \frac{(tL)^{\ell}}{\ell!} \right) \left( \sum_{\ell=0}^{p} (-1)^{\ell} a_{\ell,p} \frac{(tL)^{\ell}}{\ell!} \right)^{-1}, $$

with $a_{\ell,p} = \frac{\frac{p!}{(p-\ell)!}}{\frac{(2p)!}{(2p-\ell)!}}$. We rewrite the Padé approximant $P_{p,p}(tL)$ as

$$ \begin{aligned} P_{p,p}(tL) &= \left( P \left( \sum_{\ell=0}^{p} a_{\ell,p} \frac{(tD)^{\ell}}{\ell!} \right) P^{-1} \right) \left( P \left( \sum_{\ell=0}^{p} (-1)^{\ell} a_{\ell,p} \frac{(tD)^{\ell}}{\ell!} \right) P^{-1} \right)^{-1} \\ &= P \left( \sum_{\ell=0}^{p} a_{\ell,p} \frac{(tD)^{\ell}}{\ell!} \right) \left( \sum_{\ell=0}^{p} (-1)^{\ell} a_{\ell,p} \frac{(tD)^{\ell}}{\ell!} \right)^{-1} P^{-1} \equiv P\tilde{D}P^{-1}, \end{aligned} $$
where $\tilde{D}$ is still a diagonal matrix, as a product of two diagonal matrices. The diagonal terms of $\tilde{D}_{j,j}$ ($j = 1, \dots, d$) of $\tilde{D}$ thus writes

$$ 
\begin{aligned}
\tilde{D}_{j,j} &= \left( \sum_{\ell=0}^{p} i^{\ell} a_{\ell,p} \frac{(t\alpha_j)^{\ell}}{\ell!} \right) \left( \sum_{\ell=0}^{p} (-1)^{\ell} i^{\ell} a_{\ell,p} \frac{(t\alpha_j)^{\ell}}{\ell!} \right)^{-1} \\
&= \left( \sum_{\ell=0}^{\lfloor \frac{p}{2} \rfloor} (-1)^{\ell} a_{\ell,p} \frac{(t\alpha_j)^{2\ell}}{2\ell!} + i \sum_{\ell=0}^{\lfloor \frac{p}{2} \rfloor - 1} (-1)^{\ell} a_{\ell,p} \frac{(t\alpha_j)^{2\ell+1}}{(2\ell+1)!} \right) \\
&\quad \left( \sum_{\ell=0}^{\lfloor \frac{p}{2} \rfloor} (-1)^{\ell} a_{\ell,p} \frac{(t\alpha_j)^{2\ell}}{2\ell!} - i \sum_{\ell=0}^{\lfloor \frac{p}{2} \rfloor - 1} (-1)^{\ell} a_{\ell,p} \frac{(t\alpha_j)^{2\ell+1}}{(2\ell+1)!} \right)^{-1} = \lambda_j^{p,+} / \lambda_j^{p,-}.
\end{aligned}
$$

Observing $\lambda_j^{p,+} = \overline{\lambda_j^{p,-}}$, $\forall j = 1, \dots, d$ enables to deduce that $\text{Sp}(P_{p,p}(tL)) \subset \mathcal{C}(0, 1)$. Now let $t \in \mathbb{R}$ and $n \in \mathbb{N}$ and observe that we get similarly $(P_{p,p}(tL))^n = P \tilde{D}^n P^{-1}$, where $\tilde{D}^n$ is again a diagonal matrix with diagonal terms that belongs to the unit circle, thus a unitary matrix. Therefore $\|\tilde{D}^n\| = 1$ and the required upper bound (4.3) follows with the constant $C = \text{cond}(P) = \|P\|\|P^{-1}\|$. $\square$

We end up this section by making some remarks.

**Remark 4.2.** For the case $p \neq q$, it is possible to write down the euclidean division of $h_{p,q} \in \mathbb{R}^p[iX]$ by $k_{p,q} \in \mathbb{R}^q[iX]$. There exist $Q_{p-q} \in \mathbb{R}^{p-q}[iX]$ and $R \in \mathbb{R}^l[iX]$ with $l < q$ such that

$$ h_{p,q}(z) = k_{p,q}(z)Q_{p-q}(z) + R(z). $$

Hence, we have

$$ P_{p,q}(z) = \frac{h_{p,q}(z)}{k_{p,q}(z)} = Q_{p-q}(z) + \frac{R(z)}{k_{p,q}(z)}. $$

In the examples presented in Section 3, the argument $z$ is a pure imaginary number such that $|z| \to +\infty$ when $N_x \to +\infty$. Hence, when $p > q$, one directly obtains that $|Q_{p-q}(z)| \to +\infty$ as $|z| \to +\infty$ which leads to severe stability issues making this case usefulness in practice. A specific example is $P_{2,1}(z)$ for which we can prove that $|P_{2,1}(z)| > 1$ for all $z \in i\mathbb{R}^*$. We will present a numerical illustration in Section 7.

The case $p < q$ turns out to be more complicated. Indeed, if $|p-q|$ is large enough ($|p-q| > 2$ typically), instabilities are observed in the sense that (4.3) is not satisfied. Indeed, we investigate numerically the limit when $n \to +\infty$ of $(P_{p,q}(zT/n))^n$ for a fixed final time $T$ and with $z \in i\mathbb{R}$ such that $|z|$ is large. First, we observe numerically that even if $|P_{p,q}(z)| \to 0$ as $|z| \to +\infty$, some values of $|P_{p,q}(z)|$ are larger than 1, thus we observe that $|(P_{p,q}(zT/n))^n| \to +\infty$ as $n \to \infty$ for a fixed final time $T = 10$. Finally, still for the case $p < q$ but when $1 \le |p-q| \le 2$, we observed that the method is stable $|P_{p,q}(z)| \le 1, \forall z \in i\mathbb{R}$.

**Remark 4.3.** Another simple approximation is the Taylor expansion of the exponential, which corresponds to a Padé approximant whose degree denominator $q$ is zero:

$$ \exp(tL) \approx P_{p,0}(tL) = \sum_{\ell=0}^{p} \frac{(tL)^{\ell}}{\ell!}, $$

for which a similar error estimate can be easily derived $\|\exp(tL) - P_{p,0}(tL)\| \le Ct^{p+1}$ but, in that case, as discussed in the previous remark, stability issues are observed.
**5. Modified Lawson methods.** In this section, we perform a numerical analysis of the modified Lawson schemes, in which the exponential of the matrix $\exp(tL)$ is approximated by $\mathfrak{e}\mathfrak{r}\mathfrak{p}(tL)$ verifying for all $t \in \mathbb{R}$ and all $n \in \mathbb{N}: \|(\mathfrak{e}\mathfrak{r}\mathfrak{p}(tL))^{n}\| \le C$ (which is the case for the Padé approximant due to (4.3)). To do so, we consider the following system of ODEs (where $g$ is a smooth function which is assumed to not depend on $t$ for simplicity)

$$ (5.1) \qquad \dot{u}(t) = Lu(t) + g(u(t)), \quad u(0) = u_0 \in \mathbb{R}^d. $$

General explicit RK schemes with modified Lawson method are considered latter on, but we first give a proof of convergence for the following first order scheme

$$ (5.2) \qquad \tilde{u}^{n+1} = \mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL)\tilde{u}^n + \Delta t\mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL)g(\tilde{u}^n), \quad u^0 = u_0 \in \mathbb{R}^d. $$

We perform the following assumptions on the nonlinearity $g$ and the linear part $L$.

**ASSUMPTION 5.1. Assumptions on $g$**

* $g$ is a locally Lipschitz continuous function from $\mathbb{R}^d$ to $\mathbb{R}^d$,

* $g'$ is continuous and locally bounded from $\mathbb{R}^d$ to $\mathbb{R}^d$.

These assumptions ensure that (5.1) has a unique solution $u$ defined on some interval $[0, T]$. From the property on $g$, we deduce $\|g(u(t))\| \le C(\|u(t)\| + 1), \forall t \in [0, T]$ where the constant $C > 0$ may depend on $T$, or say on a given compact set of $\mathbb{R}^d$ that contains the range $u([0, T])$ and $\{0\}$.

**Remark 5.1.** A more tractable assumption would be to consider a global Lipschitz continuous function $g$ so that the above estimate is still available but now at any point $u \in \mathbb{R}^d$ with a uniform constant $C$: $\|g(u)\| \le C(\|u\| + 1), \forall u \in \mathbb{R}^d$. In order to lighten the next proofs, we will use that global inequality and discuss in a subsequent remark how to adapt the proofs to the local Lipschitz case only.

**ASSUMPTION 5.2. Assumptions on $L$**

* $L$ is conjugated to a skew Hermitian matrix: $L = PHP^{-1} (H^* = -H)$,

* $\exists C > 0$ such that $\mathrm{cond}(P) \equiv \|P\|\|P^{-1}\| \le C$.

In the proofs, we will use the following version of the Gronwall lemma.

**LEMMA 5.2 (Gronwall Lemma).** Let $(a_n), (b_n)$ be two nonnegative sequences and a constant $M \ge 0$ satisfying $a_{n+1} \le M + \sum_{k=0}^n a_k b_k$ for any natural number $n$, then $a_{n+1} \le M \exp(\sum_{k=0}^n b_k)$.

**5.1. First order.** In this subsection, we prove the convergence of the first order modified Lawson methods introduced above. For the first order scheme (5.2), we are able to prove the following theorem.

**THEOREM 5.3.** Under the assumptions 5.1 and 5.2 on $g$ and $L$ stated above and with $\mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL)$ such that $\|\exp(\Delta tL) - \mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL)\| \le C\Delta t^{r+1}$, the first order modified Lawson method (5.2) satisfies the following error estimate

$$ \|u(t^n) - \tilde{u}^n\| \le C\Delta t^{\min(1,r)}, \quad \forall n \le N, \text{ with } N \text{ such that } N\Delta t \le T. $$

*Proof.* The first order modified Lawson scheme writes

$$ (5.3) \qquad \tilde{u}^{n+1} = \mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL)\tilde{u}^n + \Delta t\mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL)g(\tilde{u}^n), $$

and by induction, it comes

$$ (5.4) \qquad \tilde{u}^n = (\mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL))^n \tilde{u}^0 + \Delta t \sum_{k=0}^{n-1} (\mathfrak{e}\mathfrak{r}\mathfrak{p}(\Delta tL))^{n-k} g(\tilde{u}^k). $$
One thus deduces that $\tilde{u}^n$ is bounded by considering the norm of (5.4) to get

$$ \|\tilde{u}^n\| \le C\|\tilde{u}^0\| + C\Delta t \sum_{k=0}^{n-1} (\|\tilde{u}^k\| + 1), $$

and we conclude thanks to the Gronwall lemma. In comparison with (5.3), the exact solution writes

$$ \begin{aligned} u(t^{n+1}) &= \exp(\Delta t L)u(t^n) + \Delta t \int_{t^n}^{t^{n+1}} e^{(t^{n+1}-t)L} g(u(t))dt \\ (5.5) \qquad \qquad \qquad \qquad \qquad \qquad \qquad &= \exp(\Delta t L)u(t^n) + \Delta t e^{\Delta t L} g(u(t^n)) + \mathbf{r}^n, \end{aligned} $$

where $\mathbf{r}^n$ denotes the local truncation error

$$ \begin{aligned} \mathbf{r}^n &= \int_{t^n}^{t^{n+1}} (t-t^n) \int_0^1 \frac{d}{dt} (e^{(t^{n+1}-t)L} g(u(t)))_{t=t^n+\sigma(t-t^n)} d\sigma dt \\ &= \int_{t^n}^{t^{n+1}} (t-t^n) \int_0^1 \exp(\Delta t - \sigma(t-t^n)) [-Lg(u(\cdot)) + g'(u(\cdot)) \frac{du}{dt}(\cdot)]_{t^n+\sigma(t-t^n)} d\sigma dt. \end{aligned} $$

Defining the local error as $e^{n+1} = u(t^{n+1}) - \tilde{u}^{n+1}$ and considering the difference of the two above equations (5.3) and (5.5), one obtains

$$ \begin{aligned} e^{n+1} &= \exp(\Delta t L)u(t^n) - \mathfrak{crp}(\Delta t L)\tilde{u}^n + \Delta t \left[ \exp(L\Delta t)g(u(t^n)) - \mathfrak{crp}(\Delta t L)g(\tilde{u}^n) \right] + \mathbf{r}^n \\ &= \exp(\Delta t L)e^n + \Delta t \exp(\Delta t L) [g(u(t^n)) - g(\tilde{u}^n)] \\ &\quad + (\exp(\Delta t L) - \mathfrak{crp}(\Delta t L)) [\tilde{u}^n + \Delta t g(\tilde{u}^n)] + \mathbf{r}^n. \end{aligned} $$

By induction, we deduce from there

$$ e^n = \exp(n\Delta t L)e^0 + \Delta t \sum_{k=0}^{n-1} \exp((n-k)\Delta t L)\delta_k + \sum_{k=0}^{n-1} \exp((n-1-k)\Delta t L)(\gamma_k + \mathbf{r}^k), $$

with $\delta_k = [g(u(t^k)) - g(\tilde{u}^k)]$, $\gamma_k = (\exp(\Delta t L) - \mathfrak{crp}(\Delta t L)) [\tilde{u}^n + \Delta t g(\tilde{u}^n)]$. First, the local truncation error satisfies $\|\mathbf{r}^n\| \le C\Delta t^2$. Second, using the upper bound $\|\exp(Lt)\| \le C, \forall t \in [0, T]$ and the Lipschitz bound $\mathrm{Lip}(g)$ on $g$, we deduce $\|\exp((n-k)\Delta t L)\delta_k\| \le C\mathrm{Lip}(g)\|e^k\|$. Finally, using the assumption $\|\exp(\Delta t L) - \mathfrak{crp}(\Delta t L)\| \le C\Delta t^{r+1}$ together with $\|\tilde{u}^n\| \le C, \forall n \le N$ we get

$$ \begin{aligned} \|e^n\| &\le \|\exp(n\Delta t L)\| \|e^0\| + C\mathrm{Lip}(g)\Delta t \sum_{k=0}^{n-1} \|e^k\| + Cn(\Delta t^2 + \Delta t^{r+1}) \\ &\le C\mathrm{Lip}(g)\Delta t \sum_{k=0}^{n-1} \|e^k\| + CT(\Delta t + \Delta t^r). \end{aligned} $$

We can conclude using again the Gronwall lemma 5.2.

$$ (5.6) \qquad \|e^n\| \le CT(\Delta t + \Delta t^r) \exp(C\mathrm{Lip}(g) \sum_{k=0}^{n-1} \Delta t), $$

whence finally $\|e^n\| \le C\Delta t^{\min(1,r)}$, provided that $0 \le n\Delta t \le T$. □
**Remark 5.4 (Local Lipschitz continuity).** When considering a only local Lipschitz continuous function $g$, the previous proof has to be slightly adapted, proving the boundedness of the numerical solution within its convergence. For the sake of completeness, we give hereafter a short proof in that case for the first order modified Lawson method. Actually, the idea of that proof generalizes without difficulty for high order modified Lawson methods discussed later on but then without such observations.

Let us consider a solution $t \mapsto u(t)$ defined for $t \in [0, T]$ and introduce the following compact set in $\mathbb{R}^d$: $K_T := \bigcup_{t \in [0, T]} \bar{B}(u(t), 1)$ which contains the range of the exact solution. Let us denote now $\text{Lip}(g) > 0$ the Lipschitz constant of $g$ on $K_T$. The above analysis is available to get (5.6), provided that any of the terms $(\tilde{u}_k)_{0 \le k \le n}$ belongs to $K_T$. Now to prove that, we proceed by induction on the integer $n$ and adjust a priori the value $\Delta t$ to be sufficiently small. To that aim, from the definition of $K_T$ the following inequality is sufficient:

$$2CT\Delta t \exp(C\text{Lip}(g)T) < 1.$$

**5.2. Arbitrary order and stage.** In this subsection, we extend the convergence of the modified Lawson Runge-Kutta method to the general case (order $m$ and $s$ stages). Of course, additional regularity assumptions have to be performed on $g$ to ensure the convergence of the Lawson Runge-Kutta method, i.e.

**ASSUMPTION 5.3. Assumptions on $g$**

* $g$ is a locally Lipschitz continuous function from $\mathbb{R}^d$ to $\mathbb{R}^d$,

* $g^{(p)}$, $p \le m + 1$ is continuous and locally bounded from $\mathbb{R}^d$ to $\mathbb{R}^d$.

The same assumptions as in the first order case are performed on $g$ and $L$.

Thus we have the following theorem.

**THEOREM 5.5.** Under assumptions 5.3 and 5.2 on $g$ and $L$ stated above and with $\text{erp}(\Delta tL)$ such that $\|\exp(\Delta tL) - \text{erp}(\Delta tL)\| \le C\Delta t^{r+1}$, the modified Lawson method M-LRK($s,m$) (2.4) of order $m$ with $s$ stages satisfies the following error estimate

$$\|u(t^n) - \tilde{u}^n\| \le C\Delta t^{\min(m,r)}, \quad \forall n \le N, \text{ with } N \text{ such that } N\Delta t \le T.$$

*Proof.* In a first step, we prove the boundedness of the solution $\tilde{u}^n$. For a given integer $n$, we first consider the stages iterations for $1 \le i \le s$, defined from

$$\tilde{u}_n^{(i)} = \text{erp}(c_i \Delta t L) \tilde{u}^n + \Delta t \sum_{j=1}^{i-1} a_{ij} \text{erp}((c_i - c_j) \Delta t L) g(\tilde{u}_n^{(j)}).$$

From the bound (4.3) and by iteration on $i$, we obtain the following estimate

$$\|\tilde{u}_n^{(i)}\| \le C(\|\tilde{u}^n\| + \Delta t).$$

Now let us introduce $F_n = \sum_{j=1}^s b_j \text{erp}((1 - c_j)\Delta t L)g(\tilde{u}_n^{(j)})$ so that the global iteration over $u^n$ solves as

$$\tilde{u}^n = \text{erp}(\Delta tL)^n \tilde{u}^0 + \Delta t \sum_{k=0}^{n-1} (\text{erp}(\Delta tL))^{n-1-k} F_k.$$

Thanks to the assumptions on $g$ and (5.7), we obtain the following inequality $\|F_n\| \le C(\|\tilde{u}^n\| + 1)$ and thus $\|\tilde{u}^n\| \le C(\|\tilde{u}^0\| + T + \Delta t \sum_{k=0}^{n-1} \|\tilde{u}^k\|)$. Finally from the Gronwall
301 lemma 5.2 we deduce the estimate on $\tilde{u}^n$

$$302 \qquad (5.8) \qquad \|\tilde{u}^n\| \le C(\|\tilde{u}^0\| + T) \exp(CT).$$

303 Next, we consider the error estimate and introduce $u^n$ the numerical solution 304 generated by the Lawson method RK(m, s) (order m and s stages) for which the 305 following error estimates is known to hold (see Theorem 4.7 in [28])

$$306 \qquad (5.9) \qquad \|u(t^n) - u^n\| \le C\Delta t^m.$$

Thus, from the following decomposition of the error

$$e^n = u(t^n) - \tilde{u}^n = u(t^n) - u^n + u^n - \tilde{u}^n = u(t^n) - u^n + \mathfrak{e}^n,$$

307 we have to estimate the difference $\mathfrak{e}^n = u^n - \tilde{u}^n$ between the solution of the classical 308 Lawson Runge-Kutta scheme $(u^n)_n$ and the one from the modified Lawson Runge- 309 Kutta scheme $(\tilde{u}^n)_n$. One can write

$$310 \qquad \begin{aligned} \mathfrak{e}^{n+1} &= u^{n+1} - \tilde{u}^{n+1} \\ &= \exp(\Delta t L) u^n - \mathfrak{e} \mathfrak{r} \mathfrak{p}(\Delta t L) \tilde{u}^n \\ &\quad + \Delta t \sum_{j=1}^s b_j \left[ \exp((1-c_j)\Delta t L) g(u_n^{(j)}) - \mathfrak{e} \mathfrak{r} \mathfrak{p}((1-c_j)\Delta t L) g(\tilde{u}_n^{(j)}) \right] \\ &= \exp(\Delta t L) \mathfrak{e}^n + \Delta t \sum_{j=1}^s b_j \left[ \exp((1-c_j)\Delta t L) (g(u_n^{(j)}) - g(\tilde{u}_n^{(j)})) \right] \\ &\quad + \gamma(1) \tilde{u}^n + \Delta t \sum_{j=1}^s b_j \gamma(1-c_j) g(\tilde{u}_n^{(j)}) \end{aligned} 311 312 313 314 315$$

316 where we introduced $\gamma(z) = \exp(z\Delta t L) - \mathfrak{e} \mathfrak{r} \mathfrak{p}(z\Delta t L)$. By induction, we get

$$317 \qquad \begin{aligned} \mathfrak{e}^n &= \exp(n\Delta t L) \mathfrak{e}^0 \\ &\quad + \sum_{k=0}^{n-1} \exp(\Delta t L(n-1-k)) \left( \Delta t \sum_{j=1}^s b_j \left[ \exp((1-c_j)\Delta t L) (g(u_k^{(j)}) - g(\tilde{u}_k^{(j)})) \right] \right) \\ &\quad (5.10) \qquad + \gamma(1) \tilde{u}^n + \Delta t \sum_{j=1}^s b_j \gamma(1-c_j) g(\tilde{u}_k^{(j)}) \end{aligned} 318 319$$

320 As we can see in (5.10), the error $\mathfrak{e}_k^{(j)} = u_k^{(j)} - \tilde{u}_k^{(j)}$ ($j = 1, \dots, s$) of the internal stages 321 of the Lawson and modified Lawson methods has to be estimated. From the definition 322 of the two Lawson schemes (2.3) and (2.4), it comes

$$323 \qquad \begin{aligned} \mathfrak{e}_k^{(i)} &= u_k^{(i)} - \tilde{u}_k^{(i)} \\ &= \exp(c_i \Delta t L) u^k - \mathfrak{e} \mathfrak{r} \mathfrak{p}(c_i \Delta t L) \tilde{u}^k \\ &\quad + \Delta t \sum_{j=1}^{i-1} a_{i,j} \left[ \exp((c_i - c_j) \Delta t L) g(u_k^{(j)}) - \mathfrak{e} \mathfrak{r} \mathfrak{p}((c_i - c_j) \Delta t L) g(\tilde{u}_k^{(j)}) \right] \\ &= \exp(c_i \Delta t L) \mathfrak{e}^k + \Delta t \sum_{j=1}^{i-1} a_{i,j} \left[ \exp((c_i - c_j) \Delta t L) (g(u_k^{(j)}) - g(\tilde{u}_k^{(j)})) \right] \\ &\quad + \gamma(c_i) \tilde{u}^k + \Delta t \sum_{j=1}^{i-1} a_{i,j} \gamma(c_i - c_j) g(\tilde{u}_k^{(j)}). \end{aligned} 324 325 326 327$$
As performed at the beginning of the proof, we iterate on $i$ and then we use the Lipschitz property of $g$ on the one side, and the error estimate $\|\exp(\Delta tL) - \mathfrak{e}rp(\Delta tL)\| \le C\Delta t^{r+1}$ together with the bound (5.8) on the other side. We thus deduce the following inequality for $\mathfrak{e}_k^{(i)}$, $i = 1, \dots, s$:

$$ \|\mathfrak{e}_k^{(i)}\| \le C\|\mathfrak{e}^k\| + C\Delta t \sum_{j=1}^{i-1} a_{i,j} \|\mathfrak{e}_k^{(j)}\| + C\Delta t^{r+1}, $$

from which it can be proved by induction

$$ (5.11) \qquad \|\mathfrak{e}_k^{(i)}\| \le C\|\mathfrak{e}^k\| + C\Delta t^{r+1}, \quad i = 1, \dots, s. $$

Now, we come back to (5.10) and use (5.11) and the same lines of arguments as for the internal stages enable to get

$$ \|\mathfrak{e}^n\| \le C\|\mathfrak{e}^0\| + C\Delta t \sum_{k=0}^{n-1} \left( \sum_{j=1}^{s} \|\mathfrak{e}_k^{(j)}\| + C\Delta t^{r+1} \right) \le C\Delta t \sum_{k=0}^{n-1} \|\mathfrak{e}^k\| + C\Delta t^r, $$

and we deduce from the Gronwall lemma 5.2 that $\|\mathfrak{e}^n\| \le C\Delta t^r$. Finally, we conclude using (5.9). □

Remark 5.6. As observed for the first order method, the above proof extends to the local Lipschitz case, with the same arguments.

**6. Properties for the Vlasov-Ampère equations.** Regarding the Vlasov-Ampère case, the (modified) Lawson schemes enjoys some properties that are detailed in this part.

**6.1. Preserving the charge.** An interesting property of the Lawson and modified Lawson Runge-Kutta methods applied to the Vlasov-Ampère system is the preservation of the charge, i.e. the fact that the Poisson equation $\partial_x E = \int_{\mathbb{R}} f dv - 1$ is preserved along time without resolving it. Indeed, let consider the semi-discrete system

$$ (6.1) \qquad \frac{d}{dt} \left( \begin{pmatrix} \vec{f}_k(t) \\ \hat{E}_k(t) \end{pmatrix} \right) = L \left( \begin{pmatrix} \vec{f}_k(t) \\ \hat{E}_k(t) \end{pmatrix} \right) - \left( \begin{pmatrix} \widehat{(E D \vec{f})}_k \\ 0 \end{pmatrix} \right), $$

where $(Df)_\ell$ denotes a consistent finite difference operator applied to the sequence $(f_\ell)_\ell$. An example can be $(Df)_\ell = (f_{\ell+1} - f_{\ell-1})/(2\Delta v)$ but higher order or uncentered schemes can be used.

**PROPOSITION 6.1.** *The schemes (2.3) and (2.4) (where $\mathfrak{e}rp$ is chosen as $P_{p,p}$ the Padé approximant) applied to (6.1) preserve the following discrete Poisson equation*

$$ ik \frac{2\pi}{b} \hat{E}_k^n = \Delta v \sum_{\ell=1}^{N_v} \hat{f}_{k,\ell}^n, \quad \forall n \ge 1, \text{ for } k \ne 0 \text{ and } \hat{E}_{k=0}^n = 0, $$

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

*Proof.* First, we write down the exponential of matrices for the Lawson and modified Lawson schemes (let us remark that for simplicity we will denote $k$ instead of $k2\pi/b$ with $[0, b]$ the space interval)

$$ \exp(tL) = \begin{pmatrix} \operatorname{diag}(\exp(-ikt\vec{v})) & \mathbf{0}^T \\ \frac{i\Delta v}{k}(1 - \exp(-ikt\vec{v})) & 1 \end{pmatrix}, \quad P_{p,p}(tL) = \begin{pmatrix} \operatorname{diag}(\vec{a}) & \mathbf{0}^T \\ \frac{i\Delta v}{k}(1 - \vec{a}) & 1 \end{pmatrix} $$
with $a_\ell = P_{p,p}(-iktv_\ell)$. We can observe that the structure of the matrix $P_{p,p}(tL)$ is very close to the one of $\exp(tL)$. Thus in the sequel, we will use the notation $\beta_\ell$ (the $\ell$-th component of $\vec{\beta} \in \mathbb{R}^{N_v}$) which is equal to $\exp(-iktv_\ell)$ for the Lawson case and which is equal to $P_{p,p}(-iktv_\ell)$ in the modified case. We present the proof for the first order Lawson case (that is forward Euler) but the proof can be generalized to arbitrary explicit Runge-Kutta schemes of order $m$ and stages $s$. Considering the $N_v$ first components of (6.1) gives

$$ (6.2) \qquad \hat{f}_{k,\ell}^{n+1} = \beta_{\ell} \hat{f}_{k,\ell}^{n} - \Delta t \beta_{\ell} (\widehat{E^n(Df^n)_{\ell}})_k, \quad \ell = 1, \dots, N_v, \quad \square $$

wheres the last component of (6.1) gives (using (6.2))

$$ \begin{aligned} \hat{E}_k^{n+1} &= \hat{E}_k^n + \sum_{\ell=1}^{N_v} \frac{i\Delta v}{k} (1 - \beta_\ell) \hat{f}_{k,\ell}^n - \Delta t \sum_{\ell=1}^{N_v} \frac{i\Delta v}{k} (1 - \beta_\ell) (\widehat{E^n(Df^n)_{\ell}})_k, \\ &= \hat{E}_k^n + \sum_{\ell=1}^{N_v} \frac{i\Delta v}{k} (1 - \beta_\ell) \hat{f}_{k,\ell}^n - \Delta t \frac{i\Delta v}{k} \sum_{\ell=1}^{N_v} (\widehat{E^n(Df^n)_{\ell}})_k - \sum_{\ell=1}^{N_v} \frac{i\Delta v}{k} (\hat{f}_{k,\ell}^{n+1} - \beta_\ell \hat{f}_{k,\ell}^n), \\ &= \hat{E}_k^n - \frac{i}{k} \hat{\rho}_k^{n+1} + \frac{i}{k} \hat{\rho}_k^n, \end{aligned} $$

where we denote $\hat{\rho}_k^n = \sum_{\ell=1}^{N_v} \hat{f}_{k,\ell}^n \Delta v$ and used $\sum_{\ell=1}^{N_v} (Df^n)_{\ell} = 0$ since the discrete operator $D$ is chosen to be consistent. Thus, assuming the Poisson equation is satisfied initially $ik\hat{E}_k^0 = \hat{\rho}_k^0$, this last relation implies by induction $ik\hat{E}_k^{n+1} = \hat{\rho}_k^{n+1}$ which proves the Poisson equation is propagated.

The RK case extends easily. The last stage is for $f$ (with $b_\ell^{(j)} = \exp(-i(1-c_j)\Delta tkv_\ell)$ or $P_{p,p}(-i(1-c_j)\Delta tkv_\ell)$)

$$ \hat{f}_{k,\ell}^{n+1} = b_\ell \hat{f}_{k,\ell}^{n} - \Delta t \sum_{j=1}^{s} b_j b_\ell^{(j)} (\widehat{E^{(j)}(Df^{(j)})_\ell})_k, \quad \ell = 1, \dots, N_v, $$

and for $E$, we have

$$ \hat{E}_k^{n+1} = \hat{E}_k^n + \sum_{\ell=1}^{N_v} \frac{i\Delta v}{k} (1 - b_\ell) \hat{f}_{k,\ell}^n - \Delta t \sum_{j=1}^{s} b_j \sum_{\ell=1}^{N_v} \frac{i\Delta v}{k} (1 - b_\ell^{(j)}) (\widehat{E^{(j)}(Df^{(j)})_\ell})_k, $$

which gives the Poisson equation at time $n+1$ assuming it is satisfied at time $n$.

**6.2. Preserving the total energy.** In this part, we propose a way to modify the numerical solution to ensure total energy preservation. Usually, preserving the total energy may induce the use of implicit methods which are not simple or easy to implement and in practice the total energy is quite well preserved when high order methods are used in time. The strategy is based on an orthogonal projection technique (already employed in [20, 26] to preserve moments in collisional kinetic equations) suggested in [27] and detailed in [31] for Hamiltonian PDEs. The strategy proposed in this latter reference has two steps: first an explicit high order time integrator is used to update the solution of a reformulated system for which the invariant is quadratic, and second a projection technique is employed for which the Lagrange multiplier can be explicitly obtained. As a result, the method preserves the order of the time integrator used in the first step, and is energy-preserving.
In our Vlasov context, the energy is quadratic so that we do not have to refor-
mulate the system so as to get a quadratic invariant. But, as explained before, our
methods preserve the Poisson equation so we modify the strategy proposed in [31] to
get energy and charge preservations, still retaining the high order accuracy. In the
sequel, we present the method for the Vlasov-Ampère case.

Let denote by $(\tilde{f}^{n+1}, \tilde{E}^{n+1})$ the numerical solution obtained from the high order Lawson or modified Lawson methods proposed above and denote by $\tilde{H}^{n+1} = \sum_{j,\ell} |v_\ell|^2 \tilde{f}_{j,\ell}^{n+1} \Delta x \Delta v + \sum_i \tilde{E}_j^{n+1} \Delta x$ the associated total energy. We used the index $\ell$ for the velocity grid and index $j$ for the grid in space $x_j = j \Delta x, j = 0, \dots, N_x-1, \Delta x = b/N_x$.

We will consider the following correction

$$
(6.3) \qquad f_{j,\ell}^{n+1} = (1 + \lambda_n) \tilde{f}_{j,\ell}^{n+1}, \quad E_i^{n+1} = (1 + \lambda_n) \tilde{E}_j^{n+1},
$$

where $\lambda_n$ will be given below and should be of order of the integrator used in the first
step. Thus the energy at time $n+1$ becomes

$$
\begin{align*}
H^{n+1} &= \sum_{j,\ell} |v_\ell|^2 f_{j,\ell}^{n+1} \Delta x \Delta v + \sum_j (E_j^{n+1})^2 \Delta x \\
&= \tilde{H}^{n+1} + \lambda_n \left[ \sum_{j,\ell} |v_\ell|^2 \tilde{f}_{j,\ell}^{n+1} \Delta x \Delta v + 2 \sum_j (\tilde{E}_j^{n+1})^2 \Delta x \right] + \lambda_n^2 \sum_j (\tilde{E}_j^{n+1})^2 \Delta x \\
&= \tilde{H}^{n+1} - \frac{\delta_n}{D_n} (\gamma_n + 2\alpha_n) + \frac{\delta_n^2 \alpha_n}{D_n^2},
\end{align*}
$$

where we denote $\lambda_n = -\delta_n/D_n$, $\gamma_n = \sum_{j,\ell} |v_\ell|^2 \tilde{f}_{j,\ell}^{n+1} \Delta x \Delta v$, $\alpha_n = \sum_j (\tilde{E}_j^{n+1})^2 \Delta x$ and
$\delta_n = \tilde{H}^{n+1} - H^n$. Thus, we look for $D_n \neq 0$ such that

$$
(6.4) \qquad - \frac{\delta_n}{D_n} (\gamma_n + 2\alpha_n) + \frac{\delta_n^2 \alpha_n}{D_n^2} = -\delta_n,
$$

so that $H^{n+1} = \tilde{H}^{n+1} - \delta_n = \tilde{H}^{n+1} - (\tilde{H}^{n+1} - H^n) = H^n$. Solving (6.4) leads to
$$
D_n = \left[ \gamma_n + 2\alpha_n + \sqrt{(\gamma_n + 2\alpha_n)^2 - 4\delta_n\alpha_n} \right] / 2.
$$

PROPOSITION 6.2. The schemes (2.3) and (2.4) (where exp is chosen as P<sub>p,p</sub> the Padé approximant) with the energy correction (6.3) applied to (6.1) preserve the discrete energy: $\mathcal{H}^{n+1} = \mathcal{H}^n$ with $\mathcal{H}^n = \sum_{j,\ell} |v_\ell|^2 f_{j,\ell}^n \Delta v \Delta x + \sum_j (E_j^n)^2 \Delta x$.

Proof. The goal is to ensure that $\lambda_n = \mathcal{O}(\Delta t^p)$, which means that the additional correction does not deteriorate the order of the scheme. First, since the schemes are of order $p$, we have $(\tilde{f}^{n+1}, \tilde{E}^{n+1}) = (\tilde{f}(t^{n+1}), \tilde{E}(t^{n+1})) + \mathcal{O}(\Delta t^{p+1})$ from which we deduce $(\tilde{f}^{n+1}, \tilde{E}^{n+1}) = (f^n, E^n) + \mathcal{O}(\Delta t)$ and thus $\delta_n = H(t^{n+1}) - H^n + \mathcal{O}(\Delta t^{p+1})$. But since $H(t^{n+1}) = H(0) = H^0 = H^n$, we have $\delta_n = \mathcal{O}(\Delta t^{p+1})$.

It remains to prove that $D^n$ is well defined. Let check that $(\gamma_n + 2\alpha_n)^2 - 4\delta_n\alpha_n > 0$,

$$
(\gamma_n + 2\alpha_n)^2 - 4\delta_n\alpha_n = \gamma_n^2 + 4\alpha_n^2 + 4\gamma_n\alpha_n - 4\delta_n\alpha_n = \gamma_n^2 + 4\alpha_n^2 + 4\alpha_n(\alpha_n - \delta_n).
$$

By the definition of $\alpha_n \equiv \sum_j (\tilde{E}_j^{n+1})^2 \Delta x > 0$ and using $\delta_n = \mathcal{O}(\Delta t^{p+1})$, we deduce
that there exists $\Delta t_0 > 0$ such that for $\Delta t < \Delta t_0$, we have $\alpha_n - \delta_n > 0$. Finally, from
its definition and for $\Delta t$ small enough, $D_n$ is of order $2\mathcal{H}^n$ which is of order 1. We
conclude $\lambda_n = \delta_n/D_n = \mathcal{O}(\Delta t^{p+1})$. $\square$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Order of time accuracy for the two-dimensional linear transport problem approximated by M-LRK(3,3) with the Padé approximant $P_{p,q}$ for different values of $p, q$.</figcaption></figure>

Remark 6.3. Let us remark that since the same correction is used for $f$ and $E$ and since the above method preserves the Poisson equation between $(\tilde{f}^{n+1}, \tilde{E}^{n+1})$the corrected solution $(f^{n+1}, E^{n+1})$ also enjoys the charge preservation.

Remark 6.4. This correction approach can be extended to Vlasov-Maxwell equations in multi-dimension in space and velocity.

**7. Numerical results.** In this section, several numerical results are given to illustrate the properties of the schemes introduced before. In the sequel, LRK($m, s$) will denote the (standard) Lawson Runge-Kutta method of order $m$ and with $s$ stages whereas M-LRK($m, s$) will denote the modified version presented above, i.e. where $\exp(tL)$ is approximated by $\exp(tL)$.

**7.1. Explicit case.** First, we investigate a simple problem for which we get an exact solution. The goal is to illustrate the error estimate from Theorem 5.5. Thus, we consider the following problem satisfied by $u(t, x, y) \in \mathbb{R}$

$$ (7.1) \qquad \partial_t u + a\partial_x u + b\partial_y u = 0, \quad x, y \in \mathbb{R}^2, a, b \in \mathbb{R}, $$

with the initial condition $u(0, x, y) = u_0(x, y)$. We considered the following parameters: $a = 1$, $b = 0.75$, $(x, y) \in [-2, 2]^2$ with periodic boundary conditions. A Fourier method is used for the $y$-direction whereas finite differences are used for the $x$-direction, and a Lawson RK(3,3) is used for the time discretization. The numerical parameters are chosen as follows: $N_x = N_y = 243$, and the final time is $T = 10\sigma\Delta x$, $\Delta t \in \{\sigma\Delta x/n, n = 1, \dots, 14\}$ with $\sigma = 1.44$ (which corresponds to the CFL condition for WENO5 method coupled with RK(3,3) method).

On Figure 1, we plot the order in time of the M-LRK(3,3) methods with different order for the Padé approximant $P_{p,q}$. It illustrates numerically the order of accuracy $\Delta t^{\min(m,r)}$, ($m = 3$ here) for different values of $r = p+q+1$ that have been proven in Theorem 5.5.

**7.2. Vlasov-Ampère models.** In this subsection, we focus on Vlasov-Ampère models in one dimensional spatial and velocity directions. The case considers the classical Vlasov-Ampère model (equivalent to the Vlasov-Poisson one)

We first consider a two stream instability (TSI) test for which the initial condition is

$$ f_0(x, v) = \frac{1}{\sqrt{2\pi}} v^2 e^{-v^2/2} (1 + \alpha \cos(kx)), $$
<figure><img src="image_3.png" /><figcaption>FIG. 2. TSI test: time evolution of the electric energy (left) and of the charge and total energy (right) obtained with LRK(4,4) method with $\Delta t = 0.1$, $N_x = N_v = 128$ ($v_{max} = 14$).</figcaption></figure>

with $x \in [0, 2\pi/k]$, $k = 0.5$, $v \in [-9, 9]$ and $\alpha = 10^{-3}$. In Figure 2, the time evolution of the electric energy (left) and the invariants (total energy and charge) are displayed. We used a LRK(4,4) method (fixed time step $\Delta t = 0.1$). The electric energy enjoys an exponential behavior and a nonlinear phase, as expected. The total energy has a very good behavior since it is preserved up to $10^{-6}$ and the charge is preserved up to machine accuracy. Then, in Figure 3, the same quantities are plotted (electric, total energies and charge) but here we use the energy correction explained in the previous section. We can observe that now, both charge and total energy are preserved, without affecting the time evolution in a sensible way of the electric energy. Next, we change the underlying Runge-Kutta method to use adaptive time steps to the embedded RK method DP(4,3) from [21] (still combined with the energy correction). The tolerance is tol= $10^{-4}$ and the time step is modified if the local error $L_{[p]}^{n+1} = \|f_{[4]}^{n+1} - f_{[3]}^{n+1}\|$ is larger than $10^{-9}$ (here $f_{[4]}^{n+1}$ (resp. $f_{[3]}^{n+1}$) denotes the numerical solution obtained by the 4-th order (resp. 3-rd order) method). Indeed, if the local error is smaller than $10^{-9}$, the iteration is accepted and the next time step becomes $\Delta t_{n+1} = \max(0.5, \min(1.5, (\text{tol}/L_{[p]}^{n+1})^{1/3}))\Delta t_n$, otherwise, the iteration is rejected and is performed with a smaller time step $\Delta t = \max(0.5, (\text{tol}/L_{[p]}^{n+1})^{1/3})\Delta t_n$. Let us remark that for this approach, the Poisson equation is satisfied up to machine accuracy $10^{-14}$, so as the total mass. In Figure 4 (left), we plot the evolution of the time steps together with the electric energy (log-scale). We observe that the method automatically computes large time steps (up to $\Delta t = 3$) in the linear phase and then, since in the nonlinear phase the amplitude of the electric field is larger, the time steps become smaller (mostly between 0.2 and 0.3) to respect the CFL condition induced by the nonlinear part. On the right side of Figure 4, we plot the time evolution of the charge and total energy which is preserved up to machine accuracy for large time.

We then consider the following initial condition for the (SLD test)

$$f_0(x, v) = \frac{1}{\sqrt{2\pi}} e^{-v^2/2} (1 + 0.5 \cos(kx)), \quad x \in [0, 2\pi/k], k = 0.5.$$

The same results are displayed in Figures 5 for this Landau case. Similar comments can be performed showing the robustness and versatility (change the Runge-Kutta method, fixed or adaptive time steps, with or without energy correction) of the approach. Let us remark that same comments can be done for M-LRK methods.
<figure><img src="image_4.png" /><figcaption>FIG. 3. TSI test: time evolution of the electric energy (left) and of the charge and total energy (right) obtained with LRK(4,4) method with $\Delta t = 0.1$, $N_x = N_v = 128$, with energy correction step ($v_{max} = 14$).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. TSI test: time evolution of the electric energy (left) and of the total energy (right) obtained with Lawson DP(4,3) method combined with an adaptive time step and with energy correction. $N_x = N_v = 128$ and $v_{max} = 14$.</figcaption></figure>

### 7.3. Linearized hybrid Vlasov equation 1dx-3dv.
We consider in this part the linearized hybrid Vlasov model described in Section 3 and studied in [12, 30]. In [12], Lawson method was employed to discretize this model but the linear part was split into two parts $L = L_1 + L_2$ (where $L_2$ contains the Maxwell part). However since $\exp(tL)$ can not be computed for all $k$ and $t$ even with scientific software, the $L_2$ part (which corresponds to the Maxwell equations) was put in the nonlinear part so that we considered the following form $\partial_t U = L_1 U + (L_2 U + N(U))$ on which Lawson Runge-Kutta integrators were applied since $\exp(tL_1)$ can be computed easily. In this work, we are able to consider modified Lawson integrators on the form $\partial_t U = LU + N(U)$ by using an approximation $\exp(tL)$ of $\exp(tL)$. First, we notice that the eigenvalues of $L$ are pure imaginary and the matrix $L$ can be written as $L = PHP^{-1}$ with $P$ independent from $k$ and $H$ is skew hermitian. As such, we are in the framework of Proposition 4.1 which suggests the use of Padé approximant to approximate $\exp(tL)$. This has been implemented within a scientific computing environnement which allows for automatic code generation.

We present a test to illustrate the good behavior of the approach in this framework. Following [30, 12], the initial condition is given by

$$f_h(t = 0, z, \mathbf{v}) = \frac{\rho_h}{(2\pi)^{3/2} \bar{v}_{\parallel} \bar{v}_{\perp}^2} \exp \left( -\frac{v_z^2}{2\bar{v}_{\parallel}^2} - \frac{(v_x^2 + v_y^2)}{2\bar{v}_{\perp}^2} \right),$$
<figure><img src="image_6.png" /><figcaption>FIG. 5. SLD test: time evolution of the electric energy (left) and of charge and total energy (right) obtained with Lawson DP(4,3) method combined with an adaptive time step and with energy correction. $N_x = N_v = 128$ and $v_{max} = 14$.</figcaption></figure>

with $z \in [0, 2\pi/k]$, $k = 2$, $\bar{v}_{\parallel} = 0.2$, $\bar{v}_{\perp} = 0.6$, $\rho_h = 0.2$ and $B_x(t=0, z) = \epsilon \sin(kz)$, the other unknown ($E_x, E_y, j_{c,x}, j_{c,y}, B_y$) are zero initially. The velocity domain is truncated to $\mathbf{v} = [-3.6, 3.6] \times [-3.6, 3.6] \times [-2.4, 2.4]$ and we consider $N_x = 27, N_{v_x} = N_{v_y} = 32, N_{v_z} = 41$ for the phase space discretization.

We are interested in the time history of the following energies (magnetic energy, electric energy, energy of the cold and hot particles) whose sum is preserved with time

$$ 
\begin{aligned}
\mathcal{H}_B(t) &= \frac{1}{2} \int (B_x^2(t, z) + B_y^2(t, z)) \, dz, \quad \mathcal{H}_E(t) = \frac{1}{2} \int (E_x^2(t, z) + E_y^2(t, z) + E_z^2(t, z)) \, dz, \\
\mathcal{H}_c(t) &= \frac{1}{2\Omega_{pe}^2} \int (j_{c,x}^2(t, z) + j_{c,y}^2(t, z)) \, dz, \quad \mathcal{H}_h(t) = \frac{1}{2} \iint |\mathbf{v}|^2 f_h(t, z, \mathbf{v}) \, dv \, dz.
\end{aligned}
$$

First, we consider the form $\partial_t U = L_1 U + (L_2 U + N(U))$ and compare the modified Lawson method RK(4,4) with Padé $P_{2,2}$ with the Lawson method (used in [12]). Here $\exp(tL_1)$ and $\exp(p(tL_1)) = P_{2,2}(tL_1)$ are computed, whereas $L_2$ (which corresponds to the Maxwell equations) is considered explicit through the nonlinear term. This strategy leads to a stability condition since the Maxwell equations are solved explicitly in time. We found out that, for RK(4,4), the condition $\Delta t \le (2\sqrt{2}/\pi)\Delta z \approx 0.11$ with our numerical parameters (see [12, 33]). In Figure 6, we consider the time evolution of the electric energy (with $\Delta t = 0.05 < 0.11$ to ensure stability) for both approaches. The two methods are almost indistinguishable and are able to recover the linear phase (for which the growth rate $\gamma = 0.095$ is in a very good agreement with the analytical one). Next, we compare the approach proposed in [12] with the present approach where $\exp(tL)$ is approximated by $\exp(p(tL)) = P_{2,2}(tL)$, which does not suffer from the stability condition from the Maxwell equations since they are now taken into account in the linear part ('Maxwell inside' the linear part). As mentioned before, the first approach requires a time step $\Delta t < 0.11$ to be stable thus we use a slightly larger time step $\Delta t = 0.12$. In Figure 7, we can observe that the simulation 'Maxwell outside' (ie Maxwell equations are outside the linear part as in [12]) is unstable with $\Delta t = 0.12$ whereas when using a Padé approximation of the whole $L$ matrix ('Maxwell inside' the linear part), no instability is observed. The same simulation with $\Delta t = 0.5$ leads to very similar results which shows the advantages of the current approach. In Figure 8, we compare different Padé approximations of $\exp(tL)$ ($P_{p,q}(tL)$ with $p, q = 1, 2$) with $\Delta t = 0.1$. As mentioned in Section 4, the choice $p > q$ leads to unstable results which is illustrated in Figure 8 for $P_{2,1}$ whereas the other choice
<figure><img src="image_7.png" /><figcaption>FIG. 6. Time evolution of the electric energy. Comparison between the exact exponential and Padé $P_{2,2}$ coupled with RK(4, 4) with $\Delta t = 0.05$. Maxwell outside.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Time evolution of the electric energy. Comparison between Maxwell inside and Maxwell outside with $P_{2,2}$ coupled with RK(4,4) with $\Delta t = 0.12$.</figcaption></figure>

leads to stable results. Then, we couple our strategy with adaptive time stepping using DP(4,3) time integrator with $\exp(tL) = P_{2,2}(tL)$. We plot in Figure 9 the time evolution of different energies (electric, cold kinetic and magnetic energies, defined in (7.2)) for which the theoretical instability rate is $\gamma = 0.095$ (see [30, 12]). Since energies are displayed, the instability rate is twice the theoretical rate. A very good agreement is observed and we can look at the time history of the time steps in Figure 7.2 (right) which shows that the method is able to take very large time steps (around $\Delta t \approx 3$) in the linear phase (in which the fields in the nonlinear term are of small amplitude) whereas in the nonlinear phase, the time steps are smaller to capture the nonlinear effects, in an automatic way.

**8. Conclusion.** In this work, we studied Lawson and modified Lawson methods for the time integration of Vlasov-type equations. After a suitable phase-space discretization, a set of ODEs is obtained in which a linear and nonlinear parts can be obtained which is amenable to Lawson methods. The analysis of convergence of the modified Lawson methods is carried out and leads to a rate of convergence $O(\Delta t^{\min(m,r)})$ where $m$ is the order of the underlying Runge-Kutta method and $r$ is
<figure><img src="image_9.png" /><figcaption>FIG. 8. Time evolution of the electric energy. Comparison of different Padé approximants coupled with RK(3,3) with $\Delta t = 0.1$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 9. $P_{2,2}(tL)$ coupled with $DP(3,3)$. Time evolution of the electric energy (left), cold kinetic energy (middle) and time steps (right).</figcaption></figure>

the order of approximation of $\exp(tL)$ by $e^{\tau p}(tL)$. It turns out that Padé approximant is a suitable choice regarding its spectral properties leading to stability properties of the resulting numerical scheme. Several numerical tests emphasize the good behavior of the schemes, in particular they enable to use larger time steps since they are not constraint by the stringent CFL condition coming from the linear part. Moreover, in the Vlasov-Ampère case, since the linear part is solved exactly, it leads to a scheme that preserves the Poisson equation exactly. Let us remark that the construction of high order time integrator for multi-dimensional Vlasov-Maxwell equations can be performed using different approaches like standard Runge-Kutta or (Hamiltonian) splitting techniques. However, these approaches have some drawbacks. On the one side, standard RK approaches (like [11, 2, 3]) suffer from stability condition from the linear part and charge preservation is not ensured. On the other side, even if Hamiltonian splittings benefit from structure preserving properties when phase-space mesh is fine enough (see [14, 9]), for high dimensional problems however, the use of coarse meshes destroys the good long time behavior and the construction of high order splitting leads to very costly methods (see [12]). Hence, the proposed method turns out to be a viable alternative to design high order time integrators for Vlasov equations.
552

553

**Acknowledgments.** We would like to acknowledge the assistance of volunteers in putting together this example manuscript and supplement.

554

REFERENCES

555

[1] G. A. BAKER, *Essentials of Padé approximants*, New York, Academic Press, 1975.

556

[2] J. W. BANKS AND J. A. F. HITTINGER, *A new class of non-linear, finite-volume methods for vlasov simulation*, IEEE Transactions on Plasma Science, 38 (2010), pp. 2198-2207.

557

558

[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 (2019), pp. B953-B982.

559

560

561

[4] H. BARUCQ, M. DURUFLÉ, AND M. N'DIAYE, *High-order Padé and singly diagonally runge-kutta schemes for linear odes, application to wave propagation problems*, Numerical Methods for Partial Differential Equations, (2018), pp. 760-798, https://doi.org/10.1002/num.22228.

562

563

564

[5] H. L. BERK, B. N. BREIZMAN, AND M. PEKKER, *Numerical simulation of bump-on-tail instability with source and sink*, Phys. Plasmas, 2(8) (1995), pp. 3007-3016.

565

566

[6] J. BERNIER, F. CASAS, AND N. CROUSEILLES, *Splitting methods for rotations: Application to vlasov equations*, SIAM Journal on Scientific Computing, 42 (2020), pp. A666-A697.

567

568

[7] C. K. BIRDSALL AND A. B. LANGDON, *Plasma Physics via computer simulation*, Taylor and Francis, 2005.

569

570

[8] M. CALIARI, P. KANDOLF, A. OSTERMANN, AND S. RAINER, *Comparison of software for computing the action of the matrix exponential*, BIT Numerical Mathematics, 54 (2014), pp. 113-128.

571

572

573

[9] F. CASAS, N. CROUSEILLES, E. FAOU, AND M. MEHRENBERGER, *High-order hamiltonian splitting for vlasov-poisson equations*, Numerische Mathematik, 135 (2017).

574

575

[10] C. CHENG AND G. KNORR, *The integration of the vlasov equation in configuration space*, Journal of Computational Physics, 22 (1976), pp. 330-351.

576

577

[11] Y. CHENG, I. GAMBA, F. LI, AND P. MORRISON, *Discontinuous galerkin methods for the vlasov-maxwell equations*, SIAM Journal on Numerical Analysis, 52 (2014), pp. 1017-1049.

578

579

[12] A. CRESTETTO, N. CROUSEILLES, Y. LI, AND J. MASSOT, *Comparison of high-order eulerian methods for electron hybrid model*, Journal of Computational Physics, 451 (2022).

580

581

[13] N. CROUSEILLES, L. EINKEMMER, AND E. FAOU, *Hamiltonian splitting for the vlasov-maxwell equations*, Journal of Computational Physics, 283 (2015), pp. 224-240.

582

583

[14] N. CROUSEILLES, L. EINKEMMER, AND E. FAOU, *An asymptotic preserving scheme for the relativistic vlasov-maxwell equations in the classical limit*, Computer Physics Communications, 209 (2016), pp. 13-26.

584

585

586

[15] N. CROUSEILLES, L. EINKEMMER, AND J. MASSOT, *Exponential methods for solving hyperbolic problems with application to kinetic equations*, Journal of Computational Physics, 420 (2020).

587

588

589

[16] N. CROUSEILLES, L. EINKEMMER, AND M. PRUGGER, *An exponential integrator for the drift-kinetic model*, Computer Physics Communications, 224 (2018), pp. 144-153.

590

591

[17] N. CROUSEILLES, M. MEHRENBERGER, AND E. SONNENDRÜCKER, *Conservative semi-lagrangian schemes for the vlasov equation*, Journal of computational physics, 229 (2010), pp. 1927-1953.

592

593

594

[18] N. CROUSEILLES, M. MEHRENBERGER, AND F. VECIL, *Discontinuous Galerkin semi-Lagrangian method for Vlasov-Poisson*, ESAIM: Proceedings, (2011), p. 21.

595

596

[19] N. CROUSEILLES, P. NAVARO, AND E. SONNENDRÜCKER, *Charge conserving grid based methods for the vlasov-maxwell equations*, Comptes rendus de Mécanique, 342 (2014), pp. 636-646.

597

598

[20] G. DIMARCO AND R. LOUBERE, *Towards an ultra efficient kinetic scheme. part ii: The high order case*, Journal of computational physics, 255 (2013), pp. 699-719.

599

600

[21] J. DORMAND AND P. PRINCE, *A family of embedded runge-kutta formulae*, Journal of Computational and Applied Mathematics, 6 (1980), pp. 19-26.

601

602

[22] L. EINKEMMER, *An adaptive step size controller for iterative implicit methods*, Applied Numerical Mathematics, 132 (2018), pp. 182-204.

603

604

[23] L. EINKEMMER, *A performance comparison of semi-lagrangian discontinuous galerkin and spline based vlasov solvers in four dimensions*, Journal of Computational Physics, 376 (2019), pp. 937-951.

605

606

607

[24] L. EINKEMMER, M. TOKMAN, AND J. LOFFELD, *On the performance of exponential integrators for problems in magnetohydrodynamics*, J. Comput. Phys., 330 (2017), pp. 550-565.

608

609

[25] F. FILBET AND E. SONNENDRÜCKER, *Comparison of eulerian vlasov solvers*, Computer Physics Communications, 150 (2003), pp. 247-266.

610
611

612

613

614

615

616

617

618

619

620

621

622

623

624

625

626

627

628

629

630

631

632

633

634

635

636

637

638

639

640

641

642

643

644

645

646

647

648

649

650

651

652

653

654

655

656

657

658

659

660

661

662

663

664

665

666

667

668

669

[26] J. HAACK AND I. GAMBA, *Conservative deterministic spectral boltzmann solver near the grazing collisions limit*, 28th Rarefied Gas Dynamics Conference, (2012).

[27] E. HAIRER, G. WANNER, AND C. LUBICH, *Geometric numerical integration : structure-preserving algorithms for ordinary differential equations*, Springer, Berlin New York, 2006.

[28] M. HOCHBRUCK, J. LEIBOLD, AND A. OSTERMANN, *On the convergence of lawson methods for semilinear stiff problems*, Numerische Mathematik, 145 (2020), pp. 553–580.

[29] M. HOCHBRUCK AND A. OSTERMANN, *Exponential integrators*, Acta Numerica, 19 (2010), pp. 209–286.

[30] F. HOLDERIED, S. POSSANNER, A. RATNANI, AND X. WANG, *Structure-preserving vs. standard particle-in-cell methods: The case of an electron hybrid model*, Journal of Computational Physics, 402 (2020).

[31] C. JIANG, Y. WANG, AND Y. GONG, *Explicit high-order energy-preserving methods for general hamiltonian partial differential equations*, J. Comput Appl. Math., 388 (2021), p. 113298.

[32] A. KLIMAS AND W. FARRELL, *A splitting algorithm for vlasov simulation with filamentation filtration*, Journal of Computational Physics, 110 (1994), pp. 150–163.

[33] K. KORMANN AND E. SONNENDRÜCKER, *Energy-conserving time propagation for a structure-preserving particle-in-cell vlasov-maxwell solver*, Journal of Computational Physics, 425 (2021), p. 109890.

[34] M. KRAUS, K. KORMANN, P. MORRISON, AND E. SONNENDRÜCKER, *Gempic: geometric electromagnetic particle-in-cell methods*, Journal of Plasma Physics, 83 (2017), p. 905830401, https://doi.org/10.1017/S002237781700040X.

[35] J. D. LAWSON, *Generalized runge-kutta processes for stable systems with large lipschitz constants*, SIAM Journal on Numerical Analysis, 4 (1967), pp. 372–380, https://arxiv.org/abs/https://doi.org/10.1137/0704033.

[36] J. D. LAWSON, *An order six runge-kutta process with extended region of stability*, Journal on Numerical Analysis, (1967), https://doi.org/10.1137/0704056.

[37] M. LESUR, *The Berk-Breizman Model as a Paradigm for Energetic Particle-driven Alfvén Eigenmodes*, theses, Ecole Polytechnique X, Dec. 2010.

[38] C. LI, X. ZHU, AND C. GU, *Matrix padé-type method for computing the matrix exponential*, Applied Mathematics, 2 (2011), pp. pp. 247–253.

[39] E. MADAULE, M. RESTELLI, AND E. SONNENDRÜCKER, *Energy conserving discontinuous galerkin spectral element method for the vlasov–poisson system*, Journal of Computational Physics, 279 (2014), pp. 261–288.

[40] A. MARTÍNEZ, L. BERGAMASCHI, M. CALIARI, AND M. VIANELLO, *A massively parallel exponential integrator for advection-diffusion models*, Journal of Computational and Applied Mathematics, 231 (2009), pp. 82–91.

[41] C. MOLER AND C. VAN LOAN, *Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later*, SIAM Review, 45 (2003), pp. 3–49.

[42] J.-M. QIU AND A. CHRISTLIEB, *A conservative high order semi-lagrangian weno method for the vlasov equation*, Journal of Computational Physics, 229 (2010), pp. 1130–1149, https://doi.org/10.1016/j.jcp.2009.10.016.

[43] 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 (2011), pp. 6203–6232.

[44] N. SIRCOMBE AND T. ARBER, *Valis: A split-conservative scheme for the relativistic 2d vlasov-maxwell system*, Journal of Computational Physics, 228 (2009), pp. 4773–4788.

[45] E. SONNENDRÜCKER, J. ROCHE, P. BERTRAND, AND A. GHIZZO, *The semi-lagrangian method for the numerical resolution of the vlasov equation*, Journal of Computational Physics, 149 (1999), pp. 201–220.

[46] C. TRONCI, *Hamiltonian approach to hybrid plasma models*, Journal of Physics A: Mathematical and Theoretical, 43 (2010), p. 375501.

[47] C. TRONCI, E. TASSI, E. CAMPOREALE, AND P. J. MORRISON, *Hybrid vlasov-MHD models: Hamiltonian vs. non-hamiltonian*, Plasma Physics and Controlled Fusion, 56 (2014), p. 095008.

[48] J. P. VERBONCOEUR, *Particle simulation of plasmas: review and advances*, Plasma Physics and Controlled Fusion, (2005).

[49] C. YANG AND F. FILBET, *Conservative and non-conservative methods based on hermite weighted essentially-non-oscillatory reconstruction for vlasov equations*, Journal of Computational Physics, 279 (2014), pp. 18–36.