# Semi-Lagrangian SAV method for Vlasov-Maxwell equations

Nicolas Crouseilles, Hongtao Liu, Yukun Yue

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

## To cite this version:

Nicolas Crouseilles, Hongtao Liu, Yukun Yue. Semi-Lagrangian SAV method for Vlasov-Maxwell equations. 2025. hal-05273412

HAL Id: hal-05273412

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

Preprint submitted on 23 Sep 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Semi-Lagrangian SAV method for Vlasov-Maxwell equations

Nicolas Crouseilles*

Hongtao Liu†

Yukun Yue‡

## Abstract

In this work, a new numerical method to approximate the solution of the Vlasov-Maxwell equations is presented. The method uses a phase space discretization and its main properties are: energy and charge conservation thanks to a semi-implicit treatment of the Maxwell equations, but allowing for an explicit and efficient update of the unknown. One of the main ingredients lies in the introduction of an auxiliary scalar variable inspired from the Scalar Auxiliary Variable (SAV) approach [23] together with a suitable splitting inspired from [19] which enables the use of a semi-Lagrangian method.

## 1 Introduction

The modeling of plasmas is widely performed using fluid models but when the collisionality is weak, kinetic models are preferred. Since they take into account a phase space description, kinetic models can describe out of equilibrium plasmas. One of the famous model is the Vlasov-Maxwell (VM) system which is the main goal of this work. More precisely, our goal is to describe a new numerical method for its approximation.

Many numerical strategies have been proposed in the last decades. Among them, the Particle-In-Cell method is a very popular method: it is based on the Lagrangian approximation of the unknown and the location/velocity of the particles are solution of a differential system. In spite of its excellent trade off between its cost and the reliability of the results, it is well known that these methods suffer from numerical noise which is problematic in particular in low density regions.

Eulerian methods, which use a grid of the phase space, are able to overcome this drawback by using high order approximation of the differential operators occurred in the VM system. Such approaches may suffer from strong restriction on the time step which may lead to very costly simulation since one usually wants to reach large time in plasma physics. Semi-Lagrangian methods offer a good compromise between Lagrangian and Eulerian methods since high order accuracy can be reached to capture strong gradient and fine structure of the solution still being uniformly stable by tracking back in time the characteristics equation.

An important criteria to assess the reliability of the numerical solution is the preservation of the invariants of the VM system. The VM has infinite invariants ($L^p$ norms for example) but some of them are more relevant. For instance the preservation of the total mass, the Gauss law or the total energy are very often regarded as important criteria for numerical solvers. Indeed, the non preservation of the Gauss law usually leads to unphysical results ([5, 6, 25]) and energy preservation ensures long time stability of the numerical results. But of course, the cost of the method is also important and these invariants have to be preserved independently of the mesh. These methods are usually referred as structure-preserving methods. Another important aspect is the complexity of the method. Indeed, due to the dimensionality of the Vlasov-Maxwell system, explicit methods are preferred since implicit methods usually require iterative approaches involving linear solvers which may be too costly both in terms of CPU time and memory aspects.

Recently in [19, 18], an efficient energy-preserving semi-Lagrangian method for the Vlasov-Ampère (VA) and Vlasov-Maxwell (VM) systems. The method is based on a suitable splitting introduced by

<footnote>*Univ Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042 Rennes, France & ENS Rennes (nicolas.crouseilles@inria.fr).</footnote>

<footnote>†Center for Mathematical Plasma Astrophysics, Department of Mathematics, KU Leuven, Leuven, 3001, Belgium (hongtao.liu@kuleuven.be).</footnote>

<footnote>‡Department of Mathematics, University of Wisconsin-Madison, Madison, WI 53706 (yyue24@wisc.edu).</footnote>
[11] for the Vlasov-Ampère system combined with a midpoint time integrator; moreover, it is shown in [19, 18] it is possible to use a semi-Lagrangian method in both space and velocity that preserves exactly the total energy. The approach is generalized to VM thanks to three ingredients: (i) a new step involving the free Maxwell equations and the $v \times B$ term is considered, (ii) the free Maxwell equations are solved using a midpoint method that preserves the electromagnetic energy, (iii) the $v \times B$ term is solved using a suitable exact splitting that avoids a costly multi-dimensional semi-Lagrangian step. This results in an energy-preserving semi-Lagrangian method for the full Vlasov-Maxwell equations.

However, this approach does not preserve the Gauss law. Indeed, since the Ampère part is not solved with the transport part of the Vlasov equation, this approach is a priori not compatible with a charge-preserving approach. In [12, 13, 17], time splittings of the VM system are proposed which preserves the Gauss law by construction: indeed, solving exactly the 'kinetic' part of the splitting (which consists in the Ampère equation and the free transport part of the Vlasov equation) leads to the preservation of the Gauss law. In spite of the fact they are based on a Hamiltonian splitting which should lead to very good long time behavior, the phase space discretization breaks the symplectic properties of the time integrator. Thus, preserving both Gauss law and total energy seems to be contradictory for Eulerian methods. To overcome this apparent contradiction, we introduce an additional unknown by following the Scalar Auxiliary Variable (SAV) approach [23]. Hence, it is possible to preserve the Gauss law and the total energy for the reformulated VM system using the SAV approach.

The SAV method has been introduced in [23] for dissipative nonlinear system. It consists in a quadratization of the nonlinear free energy by introducing a new scalar variable $r(t)$ so that the free energy becomes quadratic. Hence a midpoint time integrator is used and it can be proved the new quadratic energy is preserved exactly. The extension of this approach to Hamiltonian ODE system has been proposed in [3]. As far as we know, SAV method for Hamiltonian PDE system has not been considered. Here, we introduce the SAV approach to the VA and VM system. It provides a degree of freedom that enables to preserve both the Gauss law and the total energy of the VMSAV system.

Starting from the VASAV or VMSAV system, it is possible to design a suitable time splitting based on the splitting from [12, 13, 17, 19, 18] to update the numerical unknown in an explicit fashion and thanks to the semi-Lagrangian method. On one side, the 'kinetic' part of the splitting contains the free transport and the Ampère equation so that it preserves the Gauss law as in [12, 13, 17]. On the other side, the $v$ transport term under the electric field (recast using the SAV approach) can be approximated (using [19, 18]) in such a way the total energy can be proven to be preserved exactly.

The rest of the paper is organized as follows. The new method is introduced on the Vlasov-Ampère model. Then, the Vlasov-Maxwell system is considered with a specific attention on the VM1d-2dv system. Finally extensive numerical simulations are presented.

## 2 Vlasov-Ampère system

In this section, we first consider the Vlasov-Ampère (VA) system in one dimension in space and one dimension in velocity satisfied by the distribution function $f(t, x, v) \ge 0$ and the electric field $E(t, x) \in \mathbb{R}$ with the time $t \ge 0$, the space variable $x \in [0, 1]$ and the velocity variable $v \in \mathbb{R}$. The VA system can be written as

$$ \partial_t f + v \partial_x f + E \partial_v f = 0, \quad \partial_t E = -J = - \int_{\mathbb{R}} v f dv, \qquad (1) $$

with the initial condition $f(t = 0, x, v)$ and $\partial_x E(0, x) = \int_{\mathbb{R}} f(0, x, v) dv - \int_0^1 \int_{\mathbb{R}} f(0, x, v) dv dx$. Periodic boundary conditions are considered in space and vanishing boundary conditions in the velocity direction. For this system, it is well known that the Poisson equation which is satisfied initially is actually propagated with time $\partial_x E(t, x) = \int_{\mathbb{R}} f(t, x, v) dv - \int_0^1 \int_{\mathbb{R}} f(0, x, v) dv dx$ (let recall that the total mass $m(t) = \int_0^1 \int_{\mathbb{R}} f(t, x, v) dv dx = m(0)$ is preserved with time). Moreover, the total energy defined by $\mathcal{E}(t) = \mathcal{E}_{kin}(t) + \mathcal{E}_{elec}(t)$ is preserved with time $\mathcal{E}(t) = \mathcal{E}(0)$ with

$$ \mathcal{E}_{kin}(t) = \frac{1}{2} \int_{0}^{1} \int_{\mathbb{R}} v^2 f(t, x, v) dv dx \quad \text{and} \quad \mathcal{E}_{elec}(t) = \frac{1}{2} \int_{0}^{1} E^2(t, x) dx. \qquad (2) $$
## 2.1 Splitting approaches

A very popular method to solve (1) or equivalently the Vlasov-Poisson equation is the directional splitting approach [10, 15, 21, 22, 24]. Indeed, it enables to consider linear transport problem that can be solved very efficiently using the semi-Lagrangian method. However, as we shall see, according to the way the Ampère equation is combined with the directional splitting, Gauss law or total energy preservation can be lost.

### Gauss law preserving splitting

The Gauss law splitting is the following

$$ (1)_G \partial_t f + v \partial_x f = 0 \quad \text{and} \quad \partial_t E = -J, $$

$$ (2)_G \partial_t f + E \partial_v f = 0, \quad \text{and} \quad \partial_t E = 0. $$

Solving alternatively step (1)$_G$ and (2)$_G$ will give a solver for (1) that preserves the Gauss law. Indeed, as observed in [12], each step can be solved exactly in time and this enables to prove the Poisson equation is propagated with time. However, the total energy is not preserved exactly.

Let us detail the splitting substeps and investigate the resulting total energy evolution. The substeps of the splitting (1)$_G$ and (2)$_G$ are solved following

* first step (1)$_G$:

$$ f^{(1)}(x, v) = f^n(x - v\Delta t, v), \quad E^{n+1}(x) = E^n(x) - \int_{t_n}^{t_{n+1}} \int_{\mathbb{R}} v f^n(x - v(t - t_n), v) dv dt. $$

* second step (2)$_G$:

$$ f^{n+1}(x, v) = f^{(1)}(x, v - \Delta t E^{n+1}(x)). $$

Then, one can compute the kinetic energy $\mathcal{E}_{kin}^{n+1}$

$$ \begin{aligned} \mathcal{E}_{kin}^{n+1} &= \iint \frac{v^2}{2} f^{n+1}(x, v) dx dv = \iint \frac{v^2}{2} f^{(1)}(x, v - \Delta t E^{n+1}(x)) dx dv \\ &= \iint \frac{1}{2} (w + \Delta t E^{n+1}(x))^2 f^{(1)}(x, v) dx dw \\ &= \mathcal{E}_{kin}^{(1)} + \Delta t \int E^{n+1}(x) \left( J^{(1)} + \frac{\Delta t}{2} \rho^{(1)} E^{n+1}(x) \right) dx, \end{aligned} \quad (3) $$

where $(\rho^{(1)}, J^{(1)}) = \int (1, v) f^{(1)}(x, v) dv$ denotes the density and current of the distribution function after the first step, and where the kinetic and electric energies $\mathcal{E}_{kin}^{n+1}$, $\mathcal{E}_{elec}^{n+1}$ are defined in (2). Let remark that the kinetic energy is conserved by the first step since only the $x$ variable is modified so that $\mathcal{E}_{kin}^{(1)} = \mathcal{E}_{kin}^n$. Now, let multiply (2)$_G$ by $v$ and integrate wrt $v$ to get $\partial_t J = \rho^{(1)} E^{n+1}(x)$. Observing $\rho^{(1)}$ and $E^{n+1}$ are constant in time in this step, integrating on $t \in [t_n, t_{n+1/2}]$ leads to

$$ J^{n+1/2} = J^{(1)} + \frac{\Delta t}{2} \rho^{(1)} E^{n+1}(x). $$
Hence, (3) becomes (introducing $E^{n+1/2} = (E^n + E^{n+1})/2$)

$$
\begin{align*}
\mathcal{E}_{kin}^{n+1} &= \mathcal{E}_{kin}^{n} + \Delta t \int E^{n+1} J^{n+1/2} dx \\
&= \mathcal{E}_{kin}^{n} + \Delta t \int (E^{n+1} - E^{n+1/2} + E^{n+1/2}) \left( J^{n+1/2} + \frac{E^{n+1} - E^n}{\Delta t} - \frac{E^{n+1} - E^n}{\Delta t} \right) dx \\
&= \mathcal{E}_{kin}^{n} + \Delta t \int (E^{n+1} - E^{n+1/2}) J^{n+1/2} + \Delta t \int E^{n+1/2} \left( J^{n+1/2} + \frac{E^{n+1} - E^n}{\Delta t} \right) dx \\
&\quad + \frac{1}{2} \int (E^n)^2 dx - \frac{1}{2} \int (E^{n+1})^2 dx \\
&= \mathcal{E}_{kin}^{n} + \Delta t \int (E^{n+1} - E^{n+1/2}) J^{n+1/2} dx + \frac{1}{2} \int (E^n)^2 dx - \frac{1}{2} \int (E^{n+1})^2 dx \\
&\quad + \Delta t \int E^{n+1/2} \left( J^{n+1/2} - \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \int_{\mathbb{R}} v f^n(x - v(t - t_n), v) dv dt \right) dx \\
&= \mathcal{E}_{kin}^{n} + \mathcal{E}_{elec}^{n} - \mathcal{E}_{elec}^{n+1} + \mathcal{O}(\Delta t^2),
\end{align*}
$$

since on the one side $E^{n+1} - E^{n+1/2} = \mathcal{O}(\Delta t)$ by consistency in time on $E$ and on the other side $J^{n+1/2} - \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \int_{\mathbb{R}} v f^n(x - v(t - t_n), v) dv dt = \mathcal{O}(\Delta t)$ by consistency in time from the current $J$. As this splitting is a Hamiltonian splitting, a good long time behavior is expected for the total energy (see [12]). However, the phase space discretization may actually break this property since the Hamiltonian splitting is performed at the continuous level.

## Energy-preserving splitting

The energy-preserving splitting is the following and only differs from the previous splitting in the fact that the Ampère equation is solved during the second step

$$ (1)_E \partial_t f + v \partial_x f = 0 \quad \text{and} \quad \partial_t E = 0, $$

$$ (2)_E \partial_t f + E \partial_v f = 0 \quad \text{and} \quad \partial_t E = -J. $$

Solving alternatively steps (1)$_E$ and (2)$_E$ combined with a suitable time discretization will actually give a solver for (1) that preserves the total energy (see [11, 19]).

Let us recall the time approximation of (2)$_E$ that enable to preserve the total energy exactly and investigate the Gauss law preservation. The substeps of the splitting (1)$_E$ and (2)$_E$ are solved as follows

* first step (1)$_E$

$$ f^{(1)}(x, v) = f^n(x - v\Delta t, v). $$

* second step (2)$_E$

$$ f^{n+1}(x, v) = f^{(1)}(x, v - \Delta t E^{n+1/2}(x)), \quad E^{n+1}(x) = E^n(x) - \Delta t J^{n+1/2}(x). $$

Here we used the following notations $E^{n+1/2} = (E^{n+1} + E^n)/2$, $J^{n+1/2} = J^{(1)} + \frac{\Delta t}{2} \rho^{(1)} E^{n+1/2}$, and $(\rho^{(1)}(x), J^{(1)}(x)) = \int (1, v) f^{(1)}(x, v) dv$.

For this splitting, the total energy can be proved to be conserved. Indeed, considering the kinetic energy for the updated $f^{n+1}$ gives

$$
\begin{align*}
\mathcal{E}_{kin}^{n+1} &= \iint \frac{v^2}{2} f^{n+1}(x, v) dx dv = \iint \frac{v^2}{2} f^{(1)}(x, v - \Delta t E^{n+1/2}(x)) dx dv \\
&= \iint \frac{1}{2} (w + \Delta t E^{n+1/2}(x))^2 f^{(1)}(x, v) dx dv \\
&= \mathcal{E}_{kin}^{(1)} + \Delta t \int E^{n+1/2}(x) \left( J^{(1)} + \frac{\Delta t}{2} \rho^{(1)} E^{n+1/2}(x) \right) dx.
\end{align*}
$$

Let multiply (2)$_E$ by $v$ and integrate wrt $v$ to get $\partial_t J = \rho^{(1)} E^{n+1/2}$. Observing $\rho^{(1)}$ and $E^{n+1/2}$ are constant in time during this step, integrating on $t \in [t_n, t_{n+1/2}]$ leads to

$$ J^{n+1/2} = J^{(1)} + \frac{\Delta t}{2} \rho^{(1)} E^{n+1/2}. $$
Observing in addition the kinetic energy is preserved during the first step $\mathcal{E}_{kin}^{(1)} = \mathcal{E}_{kin}^{n}$, we then get

$$ 
\begin{aligned}
\mathcal{E}_{kin}^{n+1} &= \mathcal{E}_{kin}^{n} + \Delta t \int E^{n+1/2}(x) J^{n+1/2}(x) dx \\
&= \mathcal{E}_{kin}^{n} - \Delta t \int \frac{E^{n+1}(x) + E^{n}(x)}{2} \frac{E^{n+1} - E^{n}}{\Delta t} dx \\
&= \mathcal{E}_{kin}^{n} + \mathcal{E}_{elec}^{n} - \mathcal{E}_{elec}^{n+1},
\end{aligned}
$$

where the kinetic and electric energies $\mathcal{E}_{kin}^{n+1}$, $\mathcal{E}_{elec}^{n+1}$ are defined in (2).

However, for this splitting, the Gauss law is not preserved: indeed, during the step (2)$_{E}$, the density $\rho^{(1)}$ is constant in time whereas the electric field evolves with time according to the Ampère equation. Thus, $\rho$ and $E$ can not satisfy the Gauss law. More precisely, considering the space derivative of the Ampère equation $E^{n+1} = E^n - \Delta t J^{n+1/2}$ leads to (let recall the total mass definition $m(t) = \int_0^1 \int_{\mathbb{R}} f(t, x, v) dv dx = m(0)$)

$$ 
\begin{aligned}
\partial_x E^{n+1} &= \partial_x E^n - \Delta t J^{n+1/2} + \rho^n - (\rho^n - m(0)) + (\rho^{n+1} - m(0)) - \rho^{n+1} \\
&= \partial_x E^n - (\rho^n - m(0)) + (\rho^{n+1} - m(0)) - (\rho^{n+1} - \rho^n + \Delta t J^{n+1/2}) \\
&= (\rho^{n+1} - m(0)) + \mathcal{O}(\Delta t^2),
\end{aligned}
$$

where we assumed the Gauss law is true at time $t_n$ and since $\rho^{n+1} - \rho^n + \Delta t J^{n+1/2} = \mathcal{O}(\Delta t^2)$ by consistency of the continuity equation.

## 2.2 Gauss and energy-preserving splitting: SAV formulation of Vlasov-Ampère

As we saw in the previous parts, it seems difficult to get both Gauss and energy preservation using splitting methods. Hence, in order to combine the best of the two approaches, we propose to add a degree of freedom in the problem. To do so, we consider the SAV method [23] to the (1) system. It consists in the introduction of a scalar auxiliary variable $r(t)$ which satisfies $r^2(t) = \frac{1}{2} \int_0^1 E(t, x)^2 dx$. Hence, the total energy in this case becomes quadratic $\tilde{\mathcal{E}}(t) = \frac{1}{2} \int v^2 f(t, x, v) dx dv + r^2(t)$. This quadratization procedure introduced in [23] paves the way of a suitable time integrators able to preserve the Gauss law and the total energy $\tilde{\mathcal{E}}(t)$.

The first step is to write down the SAV formulation of (1) by deriving an ODE on the auxiliary variable $r$. To do so, we differentiate in time the relation $r^2(t) = \frac{1}{2} \int_0^1 E(t, x)^2 dx$ to get

$$ 2r \frac{dr}{dt} = \int E \partial_t E dx, \quad \text{so that} \quad \frac{dr}{dt} = \frac{1}{2} \int g \partial_t E dx = -\frac{1}{2} \int g J dx, $$

where we have introduced $g(t, x) = E(t, x)/r(t)$. Hence the reformulated VASAV model becomes

$$ \partial_t f + v \partial_x f + r g \partial_v f = 0, \quad (4) $$

$$ \partial_t E = -J, \quad (5) $$

$$ \frac{dr}{dt} = -\frac{1}{2} \int g J dx. \quad (6) $$

We can easily check the preservation of the total energy $\tilde{\mathcal{E}}(t)$.

**Proposition 1.** For the system (4)-(5)-(6), the quantity $\tilde{\mathcal{E}}(t)$ is preserved with time: $\frac{d}{dt} \tilde{\mathcal{E}}(t) = 0$.

*Proof.* First, let multiply (4) by $v^2/2$ and integrate in $x, v$ to get

$$ \frac{d}{dt} \int \int \frac{v^2}{2} f dxdv - r \int \int g(t, x) v f dxdv = \frac{d}{dt} \int \int \frac{v^2}{2} f dxdv - r \int g(t, x) J dx = 0. $$

Multiplying now (6) by $r$ gives $\frac{1}{2} \frac{dr^2}{dt} = -\frac{r}{2} \int g J dx$, so that, summing the two last equations directly gives

$$ \frac{d}{dt} \tilde{\mathcal{E}}(t) = \frac{d}{dt} \left[ \int \int \frac{v^2}{2} f dxdv + r^2 \right] = 0. $$

□
Now the goal is to present a time splitting of (4)-(5)-(6). As observed previously, preserving the Gauss law requires to solve the space transport together with the Ampère equation and a midpoint strategy can be employed to preserve the total energy $E(t)$. Based on these considerations, we propose the following time splitting for (4)-(5)-(6)

$$ (1)_{GE} \quad \partial_t f + v \partial_x f = 0, \quad \partial_t E = -J, \quad \text{and} \quad \frac{dr}{dt} = 0. $$

$$ (2)_{GE} \quad \partial_t f + rg\partial_v f = 0, \quad \partial_t E = 0, \quad \text{and} \quad \frac{dr}{dt} = -\frac{1}{2} \int gJdx. $$

Now, we discuss the way we solve the two steps $(1)_{GE}$ and $(2)_{GE}$ of the new splitting.

## First subsystem

For the first subsystem $(1)_{GE}$, we follow the Gauss preserving approach by solving exactly the transport in space and the Ampère equation.

$$ f^{(1)}(x, v) = f^n(x - v\Delta t, v), $$

$$ E^{(1)}(x) = E^n(x) - \int_{t_n}^{t_{n+1}} \int_{\mathbb{R}} v f^n(x - v(t - t_n), v) dv dt $$

$$ r^{(1)} = r^n. $$

During this step, the kinetic energy is constant as well as the quadratic term $r^2$, and the Gauss law is preserved. Let reformulate the update of $(f, E, r)$ from (4) as $(f, E, r)^{n+1} = \varphi_{\Delta t}^{[1]}((f, E, r)^n)$.

## Second subsystem

For the second subsystem $(2)_{GE}$, a midpoint approximation will be used to ensure energy preservation, as in the energy preserving splitting. Hence, the update of the unknown can be written as

$$ f^{n+1}(x, v) = f^{(1)}(x, v - \Delta t r^{n+1/2} g), \quad (7a) $$

$$ E^{n+1} = E^{(1)}, \quad (7b) $$

$$ r^{n+1} - r^n = -\frac{\Delta t}{2} \int g \frac{J^{n+1} + J^{(1)}}{2} dx, \quad (7c) $$

with $r^{n+1/2} = (r^{n+1} + r^n)/2$ (since $r^{(1)} = r^n$), $J^{(1)} = \int v f^{(1)} dv$ and $g = E^{(1)}/r^n$. Let observe that since $E$ is unchanged during this step, the Gauss law is thus not affected. Let us manipulate the equations (7a)-(7b)-(7c) to show that the unknown $f$ and $r$ can be computed in an explicit fashion in spite of its implicit apparent character. Indeed, one can derive a relation on the current $J$. By multiplying (7a) by $v$ and integrating wrt $v$, we get

$$ \begin{aligned} J^{n+1}(x) &= \int v f^{n+1}(x, v) dv = \int v f^{(1)}(x, v - \Delta t r^{n+1/2} g) dv \\ &= \int (w + \Delta t r^{n+1/2} g) f^{(1)}(x, v) dv \\ &= J^{(1)}(x) + \Delta t \frac{r^{n+1} + r^n}{2} g \rho^{(1)}(x), \end{aligned} $$

where $\rho^{(1)} = \int f^{(1)} dv$. Hence, substituting $J^{n+1}$ in (7c) enables to update $r^{n+1}$ as

$$ r^{n+1} = r^n - \frac{\Delta t}{2} \int g \left( J^{(1)} + \Delta t \frac{r^{n+1} + r^n}{4} g \rho^{(1)} \right) dx, $$

which gives the following expression for $r^{n+1}$

$$ r^{n+1} = \frac{1}{1 + (\Delta t^2/8) \int g^2 \rho^{(1)} dx} \left[ r^n - \frac{\Delta t}{2} \int g \left( J^{(1)} + \frac{\Delta t}{4} r^n g \rho^{(1)} \right) dx \right]. \quad (8) $$

From the knowledge of $r^{n+1}$, it is then possible to update $f$ from (7a). Let reformulate the update of $(f, E, r)$ from (7a)-(7c) and (8) as $(f, E, r)^{n+1} = \varphi_{\Delta t}^{[2]}((f, E, r)^n)$.
The first order splitting scheme thus writes

$$ (f, E, r)^{n+1} = \varphi_{\Delta t}^{[2]} \circ \varphi_{\Delta t}^{[1]}((f, E, r)^n). \qquad (9) $$

High order versions can be obtained by composition.

## 2.3 Gauss law and energy conservations

For the scheme (9) defined previously, one can prove the following result.

**Proposition 2.** The scheme (9) preserves

* the Gauss law $\partial_x E^n = \rho^n - m(0)$ for all $n \in \mathbb{N}$,

* the following total energy $\tilde{\mathcal{E}}^n = \int \int \frac{v^2}{2} f^n(x, v) dxdv + (r^n)^2$ for all $n \in \mathbb{N}$.

*Proof.* First, we observe the first step $\varphi_{\Delta t}^{[1]}$ preserves the kinetic energy $\mathcal{E}_{kin}^{(1)} = \int \int \frac{v^2}{2} f^{(1)}(x, v) dxdv = \int \int \frac{v^2}{2} f^n(x, v) dxdv = \mathcal{E}_{kin}^n$ and of course $r^2$ which is unchanged during this first step. Let then focus on the second step $\varphi_{\Delta t}^{[2]}$. To do so, we compute the kinetic energy $\mathcal{E}_{kin}^{n+1} = \int \int \frac{v^2}{2} f^{n+1}(x, v) dxdv$ at time $t^{n+1}$

$$ \begin{aligned} \mathcal{E}_{kin}^{n+1} &= \iint \frac{v^2}{2} f^{n+1}(x, v) dxdv \\ &= \iint \frac{v^2}{2} f^{(1)}(x, v - \Delta t r^{n+1/2} g) dxdv \\ &= \frac{1}{2} \iint \left( w + \Delta t r^{n+1/2} g \right)^2 f^{(1)}(x, w) dxdw \\ &= \mathcal{E}_{kin}^{(1)} + \Delta t r^{n+1/2} \iint g \left( w + \frac{\Delta t}{2} r^{n+1/2} g \right) f^{(1)}(x, w) dxdw \\ &= \mathcal{E}_{kin}^n + \Delta t r^{n+1/2} \iint g \left( \int \left( w + \frac{\Delta t}{2} r^{n+1/2} g \right) f^{(1)}(x, w) dw \right) dx \\ &= \mathcal{E}_{kin}^n + \Delta t r^{n+1/2} \iint g \left( J^{(1)}(x) + \frac{\Delta t}{2} r^{n+1/2} g \rho^{(1)}(x) \right) dx \\ &= \mathcal{E}_{kin}^n + \Delta t r^{n+1/2} \iint g J^{n+1/2}(x) dx \\ &= \mathcal{E}_{kin}^n + \Delta t \frac{r^{n+1} + r^n}{2} \iint g \frac{J^{n+1} + J^n}{2} dx. \end{aligned} \qquad (10) $$

Now multiplying (7c) by $(r^{n+1} + r^n)$ gives

$$ \frac{(r^{n+1})^2 - (r^n)^2}{\Delta t} = - \frac{r^{n+1} + r^n}{2} \iint g \frac{J^{n+1} + J^n}{2} dx, $$

which gives the discrete energy preservation by adding the last equality with (10). $\square$

**Remark**

These calculations can be extended to the multi-dimensional version of Vlasov-Ampère. Indeed, it suffices to consider one additional scalar variable per electric field component, ie $r_x, r_y, r_z$ for $E_x, E_y, E_z$. The splitting can then be obtained by extending directly the one dimensional case presented above and the energy preservation can be deduced exactly from the proof of Prop. 2. $\square$

## 3 Vlasov-Maxwell 1dx-2dv (VM1d-2v)

In this section, the goal is to extend the strategy presented for the Vlasov-Ampère case to the Vlasov-Maxwell case. In this section, we consider the reduced Vlasov-Maxwell system in one space dimension and two velocity dimensions. Hence, let denote $(f(t, x, v), E(t, x) = (E_x, E_y)(t, x), B(t, x)) \in \mathbb{R}_+ \times \mathbb{R}^2 \times \mathbb{R}$
with $t \ge 0$, $x \in [0, L] \subset \mathbb{R}$ and $v \in \mathbb{R}^2$ the solution to the following Vlasov-Maxwell system

$$ \partial_t f + v_x \partial_x f + (E + B\Omega v) \cdot \nabla_v f = 0, \qquad (11a) $$

$$ \partial_t E_x = -J_x, \quad \text{with} \quad J_x = - \int_{\mathbb{R}^2} v_x f dv, \qquad (11b) $$

$$ \partial_t E_y = -J_y - \partial_x B, \quad \text{with} \quad J_y = - \int_{\mathbb{R}^2} v_y f dv, \qquad (11c) $$

$$ \partial_t B = -\partial_x E_y, \qquad (11d) $$

with $\Omega$ the 2x2 symplectic matrix. As in the VA case, preserving both the Gauss law and the total energy is antagonist. In the following, we thus introduce some additional variables in order to get a splitting scheme that preserves the Gauss law and a total energy.

## 3.1 SAV formulation of the VM1d-2v system

To extend the SAV method presented in the previous section to the VM case, we introduce two auxiliary variables: the first one $r_x$ is associated to $E_x$ and the second one $r_y$ is associated to $E_y$ such that

$$ r_{x,y}(t) = \frac{1}{2} \int E_{x,y}^2(t,x) dx \quad \text{and} \quad r_{x,y}(t) g_{x,y}(t,x) = E_{x,y}(t,x). $$

Hence the VM system (11) is reformulated into

$$ \partial_t f + v_x \partial_x f + (r_x g_x + B v_y) \partial_{v_x} f + (r_y g_y - B v_x) \partial_{v_y} f = 0, \qquad (12a) $$

$$ \partial_t E_x = -J_x, \qquad (12b) $$

$$ \partial_t E_y = -J_y - \partial_x B, \qquad (12c) $$

$$ \partial_t B = -\partial_x(r_y g_y), \qquad (12d) $$

$$ \dot{r}_x = -\frac{1}{2} \int g_x J_x dx, \qquad (12e) $$

$$ \dot{r}_y = -\frac{1}{2} \int g_y (J_y + \partial_x B) dx. \qquad (12f) $$

The goal is to approximate numerically the system (12) instead of (11). The strategy lies in the fact the (scalar) additional variables provide a degree of freedom in the construction of the method so that the method will combine many interesting features (no CFL, preservation of a total energy and of the Gauss law). In the sequel, we will use the following notation for the unknown $U(t) = (f, E_x, E_y, B, r_x, r_y)(t)$.

First, let prove the total energy $\tilde{\mathcal{E}}(t)$ defined by

$$ \tilde{\mathcal{E}}(t) = \iint \frac{|v|^2}{2} f(t, x, v) dxdv + r_x^2(t) + r_y^2(t) + \frac{1}{2} \int B(t, x)^2 dx, \qquad (13) $$

is preserved by the system (12a)-(12f).

**Proposition 3.** For the system (12), the quantity $\tilde{\mathcal{E}}(t)$ defined by (13) is preserved with time: $\frac{d}{dt} \tilde{\mathcal{E}}(t) = 0$.

*Proof.* We can easily check the energy preservation by multiplying (12a) by $|v|^2/2$ and integrating in $x, v$ to get

$$ \frac{d}{dt} \int \int \frac{|v|^2}{2} f dxdv - r_x \int g_x J_x dx - r_y \int g_y J_y dx = 0, $$

Now multiplying (12e) by $r_x$ and (12f) by $r_x$ leads to

$$ \frac{d}{dt}(r_x^2 + r_y^2) = - \int r_x g_x J_x dx - \int r_y g_y J_y dx - \int r_y g_y \partial_x B dx. $$

Summing the two last equation leads to

$$ \frac{d}{dt} \int \int \frac{|v|^2}{2} f dxdv + \frac{d}{dt}(r_x^2 + r_y^2) + \int r_y g_y \partial_x B dx = 0. \qquad (14) $$
Now, let focus on the last term $\int r_y g_y \partial_x B dx$. First, we integrate it by part and use the Faraday equation (12d) to get

$$\int r_y g_y \partial_x B dx = - \int \partial_x (r_y g_y) B dx = \int (\partial_t B) B dx = \frac{1}{2} \frac{d}{dt} \int B^2 dx.$$

Substituting this relation in (14) provides the preservation of $\tilde{\mathcal{E}}(t)$. □

## 3.2 Time splitting

In this part, our goal is to propose a suitable time splitting for (12) for which the preservation of the Gauss law and the preservation of $\tilde{\mathcal{E}}$ defined by (13) can be proved. The splitting will be performed into three steps that are detailed below. Let recall the notation $U(t) = (f, E_x, E_y, B, r_x, r_y)(t)$ and its discrete counterpart $U^n$.

### 3.2.1 First subsystem

To ensure Gauss law preservation, the Ampère equation has to be solved exactly in combination with the space transport part. Hence we get

$$\partial_t f + v_x \partial_x f = 0, \qquad (15a)$$

$$\partial_t E = -J, \qquad (15b)$$

$$\partial_t B = 0, \quad \dot{r}_x = \dot{r}_y = 0. \qquad (15c)$$

This subsystem can be solved exactly as in the VA case, see (4). Then, for $U^n = (f, E_x, E_y, B, r_x, r_y)^n$ given, we denote the solution of this first subsystem as

$$U^{n+1} = \varphi_{\Delta t}^{[1]}(U^n). \qquad (16)$$

Let recall that if the Poisson equation is satisfied at time $n$, then it will be satisfied at time $n+1$ so that the Gauss law is preserved during this step.

### 3.2.2 Second subsystem

The second subsystem is the following

$$\partial_t f + r_x g_x \partial_{v_x} f + r_y g_y \partial_{v_y} f = 0,$$

$$\partial_t B = -\partial_x(r_y g_y),$$

$$\dot{r}_x = -\frac{1}{2} \int g_x J_x dx,$$

$$\dot{r}_y = -\frac{1}{2} \int g_y (J_y + \partial_x B) dx,$$

$$\partial_t E = 0.$$

To preserve the total energy, we extend what has been done in the Vlasov-Ampère case and suggest the following midpoint semi-implicit numerical scheme

$$\partial_t f + r_x^{n+1/2} g_x \partial_{v_x} f + r_y^{n+1/2} g_y \partial_{v_y} f = 0, \qquad (18a)$$

$$B^{n+1} = B^n - \Delta t \partial_x (r_y^{n+1/2} g_y), \qquad (18b)$$

$$r_x^{n+1} = r_x^n - \frac{\Delta t}{2} \int g_x J_x^{n+1/2} dx, \qquad (18c)$$

$$r_y^{n+1} = r_x^n - \frac{\Delta t}{2} \int g_y (J_y^{n+1/2} + \partial_x B^{n+1/2}) dx, \qquad (18d)$$

$$E^{n+1} = E^n,$$

where we used $r_x^{n+1/2} = (r_x^n + r_y^n)/2$ (the same notation is used for the other unknown). Regarding (18a), a semi-Lagrangian method will be used

$$f^{n+1}(x, v) = f^n(x, v_x - \Delta t r_x^{n+1/2} g_x, v_y - \Delta t r_y^{n+1/2} g_y). \qquad (19)$$
Even if the time advance seems to be implicit, we will see it is possible to update the unknown with an explicit cost. Multiplying (19) by $v_x$ and then by $v_y$ and integrating wrt $v = (v_x, v_y)$ leads to

$$ 
\begin{aligned}
J_x^{n+1} &= \int_{\mathbb{R}^2} v_x f^{n+1}(x, v) dv = \int_{\mathbb{R}^2} v_x f^n(x, v_x - \Delta t r_x^{n+1/2} g_x, v_y - \Delta t r_y^{n+1/2} g_y) dv_x dv_y \\
&= \int_{\mathbb{R}^2} (w_x + \Delta t r_x^{n+1/2} g_x) f^n(x, w_x, w_y) dw_x dw_y \\
&= J_x^n + \Delta t r_x^{n+1/2} g_x \rho^n,
\end{aligned} 
\quad (20) $$

$$ J_y^{n+1} = \int_{\mathbb{R}^2} v_y f^{n+1}(x, v) dv = J_y^n + \Delta t r_y^{n+1/2} g_y \rho^n. \quad (21) $$

Similarly, from the Faraday equation (18b), one has

$$ B^{n+1/2} = \frac{1}{2}(B^{n+1} + B^n) = \frac{1}{2}(2B^n - \Delta t \partial_x (r_y^{n+1/2} g_y)) = B^n - \frac{\Delta t}{2} \partial_x (r_y^{n+1/2} g_y). \quad (22) $$

Substituting (20) in (18c) leads to

$$ 
\begin{aligned}
r_x^{n+1} &= r_x^n - \frac{\Delta t}{2} \int g_x \left[ J_x^n + \frac{\Delta t}{2} r_x^{n+1/2} g_x \rho^n \right] dx \\
&= r_x^n - \frac{\Delta t}{2} \int g_x J_x^n dx - \frac{\Delta t^2}{8} (r_x^{n+1} + r_x^n) \int g_x^2 \rho^n dx.
\end{aligned} 
$$

Now, substituting (21) in (18d) leads to

$$ 
\begin{aligned}
r_y^{n+1} &= r_y^n - \frac{\Delta t}{2} \int g_y \left[ J_y^n + \frac{\Delta t}{2} r_y^{n+1/2} g_y \rho^n + \partial_x B^{n+1/2} \right] dx, \\
&= r_y^n - \frac{\Delta t}{2} \int g_y J_y^n dx - \frac{\Delta t^2}{8} (r_y^{n+1} + r_y^n) \int g_y^2 \rho^n dx + \frac{\Delta t}{2} \int (\partial_x g_y) B^{n+1/2} dx,
\end{aligned} 
$$

and using (22) for $B^{n+1/2}$ in the last term enables to get

$$ 
\begin{aligned}
r_y^{n+1} &= r_y^n - \frac{\Delta t}{2} \int g_y J_y^n dx - \frac{\Delta t^2}{8} (r_y^{n+1} + r_y^n) \int g_y^2 \rho^n dx \\
&\quad + \frac{\Delta t}{2} \int (\partial_x g_y) \left(B^n - \frac{\Delta t}{2} \partial_x (r_y^{n+1/2} g_y)\right) dx, \\
&= r_y^n - \frac{\Delta t}{2} \int g_y J_y^n dx - \frac{\Delta t^2}{8} (r_y^{n+1} + r_y^n) \int g_y^2 \rho^n dx \\
&\quad + \frac{\Delta t}{2} \int (\partial_x g_y) B^n dx - \frac{\Delta t^2}{8} (r_y^{n+1} + r_y^n) \int (\partial_x g_y)^2 dx.
\end{aligned} 
$$

Then, thanks to these calculations, $r_{x,y}^{n+1}$ can be computed explicitly through

$$ 
\begin{aligned}
r_x^{n+1} &= \frac{1}{1 + \frac{\Delta t^2}{8} \int g_x^2 \rho^n dx} \left[ r_x^n - \frac{\Delta t}{2} \int g_x J_x^n dx - \frac{\Delta t^2}{8} r_x^n \int g_x^2 \rho^n dx \right], \\
r_y^{n+1} &= \frac{1}{1 + \frac{\Delta t^2}{8} \left( \int g_y^2 \rho^n dx + \int (\partial_x g_y)^2 dx \right)} \left[ r_y^n - \frac{\Delta t}{2} \int g_y J_y^n dx - \frac{\Delta t^2}{8} r_y^n \int g_y^2 \rho^n dx \right. \\
&\quad \left. + \frac{\Delta t}{2} \int (\partial_x g_y) B^n dx - \frac{\Delta t^2}{8} r_y^n \int (\partial_x g_y)^2 dx \right].
\end{aligned} 
\quad (23) $$

In this substep, once $r_{x,y}^{n+1}$ are updated explicitly from the relations (23)-(23) so that $B^{n+1}$ can be computed with (18b) and $f^{n+1}$ can be updated from (19). Hence, for $U^n = (f, E_x, E_y, B, r_x, r_y)^n$ given, we denote the numerical solution of this second subsystem by

$$ U^{n+1} = \varphi_{\Delta t}^{[2]}(U^n). \quad (24) $$
### 3.2.3 Third subsystem

The last subsystem we have to solve is the following

$$ 
\begin{aligned}
\partial_t f + B\Omega v \cdot \nabla_v f &= 0, \\
\partial_t E_y &= -\partial_x B, \\
\partial_t B &= 0, \\
\dot{r}_x &= \dot{r}_y = 0,
\end{aligned}
 $$

where we recall $\Omega$ is the 2x2 symplectic matrix

$$ \Omega = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}. $$

This subsystem can also be solved exactly in time

$$ f^{n+1}(x, v) = f^n(x, e^{\Omega \Delta t B^n} v), \qquad (26a) $$

$$ E_y^{n+1} = E_y^n - \Delta t \partial_x B^n, \qquad (26b) $$

$$ B^{n+1} = B^n, $$

$$ r_{x,y}^{n+1} = r_{x,y}^{n}. $$

Let us remark that the update of $f$ requires a two-dimensional interpolation in the $v$-space which may be costly to solve from a numerical point of view. However, as observed in [1, 2], this rotation in the velocity-space can be decomposed *exactly* into one dimensional shears, ie one dimensional advections. This decomposition follows from the following relation

$$ e^{\Omega \Delta t B^n} = \begin{pmatrix} \cos B^n \Delta t & \sin B^n \Delta t \\ -\sin B^n \Delta t & \cos B^n \Delta t \end{pmatrix} = \begin{pmatrix} 1 & \tan B^n \Delta t/2 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ -\sin B^n \Delta t & 1 \end{pmatrix} \begin{pmatrix} 1 & \tan B^n \Delta t/2 \\ 0 & 1 \end{pmatrix}, $$

from which we deduce one-dimensional advections

$$ f^{n+1}(v) = f^n(e^{\Omega \Delta t B^n} v) = S_{v_x \tan B^n \Delta t/2}^{v_y} S_{-v_y \sin B^n \Delta t}^{v_x} S_{v_x \tan B^n \Delta t/2}^{v_y} f^n(v), $$

where $S_{\alpha}^{v_x} f(v) = f^n(v_x + \alpha, v_y)$ and $S_{\alpha}^{v_y} f(v) = f^n(v_x, v_y + \alpha)$. Hence, one-dimensional semi-Lagrangian method can be used for each shear steps.

In this third step, $f$ and $E_y$ are update from (26a) and (26b). Hence, for $U^n = (f, E_x, E_y, B, r_x, r_y)^n$ given, we denote the solution of this third subsystem by

$$ U^{n+1} = \varphi_{\Delta t}^{[3]}(U^n). \qquad (27) $$

#### Remark

Here, we followed the splitting from [12] but other splitting can be used as the one proposed in citelapenta2 in which the free Maxwell equations are solved with the magnetic transport (26a). $\square$

### 3.3 Splitting scheme

For $U^n = (f, E_x, E_y, B, r_x, r_y)^n$ given at iteration $n$, the solution of the numerical scheme presented in the previous subsection can be written as

$$ U^{n+1} = \varphi_{\Delta t}^{[3]} \circ \varphi_{\Delta t}^{[2]} \circ \varphi_{\Delta t}^{[1]}(U^n), \qquad (28) $$

where $\varphi_{\Delta t}^{[i], i = 1, 2, 3}$ are given by (16), (24) and (27). The scheme (28) is referred as a Lie splitting, which is first order accurate in time. High order scheme can be obtained by composition (see [1]); as a typical example, the Strang splitting (which is second order accurate) would be

$$ U^{n+1} = \varphi_{\Delta t/2}^{[3]} \circ \varphi_{\Delta t/2}^{[2]} \circ \varphi_{\Delta t}^{[1]} \circ \varphi_{\Delta t/2}^{[2]} \circ \varphi_{\Delta t/2}^{[3]}(U^n). $$
## 3.4 Properties of the new splitting

One can prove the following result

**Proposition 4.** The scheme (28) preserves

* the Gauss law $\partial_x E_x^n = \rho^n - m(0)$ for all $n \in \mathbb{N}$,

* the following total energy $\tilde{\mathcal{E}}^n = \int \int \frac{|v|^2}{2} f^n(x, v) dxdv + (r_x^n)^2 + (r_y^n)^2 + \frac{1}{2} \int (B^n)^2 dx$ for all $n \in \mathbb{N}$.

*Proof.* The preservation of the Gauss law directly follows the exact integration of the first step (16) since $E_x$ is unchanged in the second and third steps.

Then, one wants to verify the energy conservation of the numerical scheme (28). Let denote $\mathcal{E}_{x,y}^n = \int_{\mathbb{R}^2 \times [0,L]} \frac{v_x^2}{2} f^n(x, v) dxdv$ and $\mathcal{E}_{kin}^n = \mathcal{E}_x^n + \mathcal{E}_y^n = \int_{\mathbb{R}^2 \times [0,L]} \frac{|v|^2}{2} f^n(x, v) dxdv$.

$$ \begin{aligned} \mathcal{E}_{kin}^{n+1} &= \mathcal{E}_x^{n+1} + \mathcal{E}_y^{n+1} \\ &= \iint \frac{v_x^2}{2} f^{n+1}(x, v) dxdv + \iint \frac{v_y^2}{2} f^{n+1}(x, v) dxdv. \end{aligned} $$

Since the first and third steps $\varphi^{[1,3]}$ given by (16) and (27) preserve the kinetic energy, it remains to study the second step $\varphi^{[2]}$ given by (18). By a small abuse of notation, this second step will consider the initial condition at iteration $n$ to update the solution at iteration $n+1$.

Multiplying (18a) by $v_x^2/2$ and integrating on $v$ give for the first term $\int \int \frac{v_x^2}{2} f^{n+1}(x, v) dxdv$

$$ \mathcal{E}_x^{n+1} = \mathcal{E}_x^n + \frac{\Delta t}{4} (r_x^{n+1} + r_x^n) \int g_x(J_x^{n+1} + J_x^n)dx, \quad (29) $$

whereas for the second term $\mathcal{E}_y^{n+1}$, we multiply (18a) by $v_y^2/2$ and integrate on $v$ to get

$$ \begin{aligned} \mathcal{E}_y^{n+1} &= \iint \frac{v_y^2}{2} f^n(x, v_x - \Delta t r_x^{n+1/2} g_x, v_y - \Delta t r_y^{n+1/2} g_y) dxdv \\ &= \frac{1}{2} \iint (w_y + \Delta t r_y^{n+1/2} g_y)^2 f^n(x, w) dxdw \\ &= \mathcal{E}_y^n + \iint \Delta t r_y^{n+1/2} g_y \left(w_y + \frac{\Delta t^2}{2} r_y^{n+1/2} g_y\right) f^n(x, w) dxdw \\ &= \mathcal{E}_y^n + \Delta t r_y^{n+1/2} \int g_y J_y^{n+1/2} dx, \end{aligned} \quad (30) $$

since $J_y^{n+1/2} = J_y^n + \frac{\Delta t}{2} r_y^{n+1/2} g_y \rho^n$. But multiplying (18c) by $2r_x^{n+1/2} = (r_x^{n+1} + r_x^n)$ gives

$$ (r_x^{n+1})^2 - (r_x^n)^2 = -r_x^{n+1/2} \Delta t \int g_x J_x^{n+1/2} dx, $$

and (18d) by $2r_y^{n+1/2} = (r_y^{n+1} + r_y^n)$ enables to get

$$ (r_y^{n+1})^2 - (r_y^n)^2 = -r_y^{n+1/2} \Delta t \int g_y J_y^{n+1/2} - r_y^{n+1/2} \Delta t \int g_y \partial_x B^{n+1/2} dx. $$

Adding now the two last equalities with (29) and (30) enables to get the following identity

$$ (r_x^{n+1})^2 + (r_y^{n+1})^2 + \mathcal{E}_x^{n+1} + \mathcal{E}_y^{n+1} = (r_x^n)^2 + (r_y^n)^2 + \mathcal{E}_x^n + \mathcal{E}_y^n - r_y^{n+1/2} \Delta t \int g_y \partial_x B^{n+1/2} dx. \quad (31) $$

It remains to consider the last term to make appear the magnetic energy. As in the continuous case, we perform an integration by parts and use the update of $B$ from the Faraday equation (18b)

$$ \begin{aligned} -r_y^{n+1/2} \Delta t \int g_y \partial_x B^{n+1/2} dx &= \Delta t \int \partial_x (r_y^{n+1/2} g_y) B^{n+1/2} dx \\ &= \frac{1}{2} \int (B^n - B^{n+1})(B^n + B^{n+1}) dx \\ &= \frac{1}{2} \int (B^n)^2 dx - \frac{1}{2} \int (B^{n+1})^2 dx. \end{aligned} \quad (32) $$
Hence one can conclude by substituting (32) in (31) to get $\tilde{\mathcal{E}}^{n+1} = \tilde{\mathcal{E}}^n$.

Hence, the total energy $\tilde{\mathcal{E}}^n$ is preserved along the second step. As mentioned before, the kinetic energy is preserved along the first and third step, whereas $B$ and $r_x, r_y$ are not updated through these first and second step, so we conclude on the preservation of $\tilde{\mathcal{E}}^n$ by the scheme (28). $\square$

# 4 Vlasov-Maxwell system

We now consider the general full Vlasov-Maxwell case satisfied by $(f(t, x, v), E(t, x), B(t, x)) \in \mathbb{R}_+ \times \mathbb{R}^3 \times \mathbb{R}^3$ with $x = (x, y, z) \in \mathbb{T} \subset \mathbb{R}^3$ and $v \in \mathbb{R}^3$. The Vlasov-Maxwell system writes as

$$ \partial_t f + v \cdot \nabla_x f + (E + v \times B) \cdot \nabla_v f = 0, \qquad (33a) $$

$$ \partial_t E = \nabla_x \times B - J, \quad \text{with} \quad J = -\int_{\mathbb{R}^3} v f dv, \qquad (33b) $$

$$ \partial_t B = -\nabla_x \times E. \qquad (33c) $$

It is well known the system (33)- preserves the total energy $\mathcal{E}(t)$ which can be written as $\frac{d}{dt}\mathcal{E}(t) = 0$ with $\mathcal{E}(t)$ given by

$$ \mathcal{E}(t) = \iint \frac{|v|^2}{2} f dx dv + \frac{1}{2} \int |E|^2 dx + \frac{1}{2} \int |B|^2 dx. $$

## 4.1 SAV formulation of the VM system

We follow the same presentation as in the VM1d-2v case. Indeed, to derive the SAV formulation of the VM system (33), we introduce the following auxiliary variables

$$ r_{x,y,z}(t)^2 = \frac{1}{2} \int E_{x,y,z}^2(t,x) dx \quad \text{and} \quad r_{x,y,z}(t) g_{x,y,z}(t,x) = E_{x,y,z}(t,x). $$

Hence the VM system (33a)-(33b)-(33c) can be reformulated into the following SAV-VM system

$$ \left\{ \begin{array}{l} \partial_t f + v \cdot \nabla_x f + (r_x g_x + v_y B_z - v_z B_y) \partial_{v_x} f \\ \qquad + (r_y g_y - v_x B_z + v_z B_x) \partial_{v_y} f + (r_z g_z + v_x B_y - v_y B_x) \partial_{v_z} f = 0, \end{array} \right. \qquad (34a) $$

$$ \partial_t E = \nabla_x \times B - J, \qquad (34b) $$

$$ \partial_t B = -\nabla_x \times u, \quad \text{with } u = (r_x g_x, r_y g_y, r_z g_z) \in \mathbb{R}^3 \qquad (34c) $$

$$ \dot{r}_x = \frac{1}{2} \int g_x(-J_x + \partial_y B_z - \partial_z B_y) dx, \qquad (34d) $$

$$ \dot{r}_y = \frac{1}{2} \int g_y(-J_y - \partial_x B_z + \partial_z B_x) dx, \qquad (34e) $$

$$ \dot{r}_z = \frac{1}{2} \int g_z(-J_z + \partial_x B_y - \partial_y B_x) dx. \qquad (34f) $$

To clarify the notation used in (34c), let write the equation for each component

$$ \partial_t B_x = -(\partial_y(r_z g_z) - \partial_z(r_y g_y)), $$

$$ \partial_t B_y = (\partial_x(r_z g_z) - \partial_z(r_x g_x)), $$

$$ \partial_t B_z = -(\partial_x(r_y g_y) - \partial_y(r_x g_x)). $$

As done in the previous section, let prove the total energy $\tilde{\mathcal{E}}(t)$ defined by

$$ \tilde{\mathcal{E}}(t) = \iint \frac{|v|^2}{2} f dx dv + \frac{1}{2} \int |B|^2 dx + (r_x^2 + r_y^2 + r_z^2), \qquad (35) $$

is preserved by the system (34a)-(34f).

**Proposition 5.** For the system (34a)-(34f), the quantity $\tilde{\mathcal{E}}(t)$ defined by (35) is preserved with time: $\frac{d}{dt}\tilde{\mathcal{E}}(t) = 0$.
Proof. To prove the energy conservation, we multiply (34a) by $|v|^2/2$ and integrate wrt x, v to get (re-
calling the notation $u = (r_x g_x, r_y g_y, r_z g_z)$)

$$
\frac{d}{dt} \int \int \frac{|v|^2}{2} f dx dv - \int u \cdot J dx = 0. \qquad (36)
$$

Multiplying (34d) by $r_x$, (34e) by $r_y$ and (34f) by $r_z$ provide (still using the notation $u_* = r_*g_*$)

$$
\begin{align}
\frac{d}{dt}(r_x^2 + r_y^2 + r_z^2) &= \int r_x g_x (-J_x + \partial_y B_z - \partial_z B_y) dx + \int r_y g_y (-J_y - \partial_x B_z + \partial_z B_x) dx \nonumber \\
&\quad + \int r_z g_z (-J_z + \partial_x B_y - \partial_y B_x) dx, \nonumber \\
&= - \int u \cdot J dx + \int u \cdot (\nabla_x \times B) dx. \tag{37}
\end{align}
$$

Summing (36) and (37) leads to

$$
\frac{d}{dt} \int \int \frac{|v|^2}{2} f dx dv + \frac{d}{dt} (r_x^2 + r_y^2 + r_z^2) = \int u \cdot (\nabla_x \times B) dx.
$$

Integrating by parts the right hand side gives $\int u \cdot (\nabla_x \times B)dx = \int (\nabla_x \times u) \cdot Bdx$ and using (34c) gives

$$
\frac{d}{dt} \int \int \frac{|v|^2}{2} f dx dv + \frac{d}{dt} (r_x^2 + r_y^2 + r_z^2) = - \int (\partial_t B) \cdot B dx = - \frac{1}{2} \frac{d}{dt} \int |B|^2 dx, \quad (38)
$$

which ends the proof.

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

4.2 Time splitting

The goal of this subsection is to present a time splitting for the VM-SAV model (34a)-(34f), for which
the Gauss law $\nabla_x \cdot E - (\rho - m(0))$ and the total energy $\tilde{\mathcal{E}}$ given by (35) are preserved.

4.2.1 First subsystem

Inspired by the splittings presented before, the first subsystem is

$$
\begin{align*}
\partial_t f + v \cdot \nabla_x f &= 0, \\
\partial_t E &= -J, \\
\partial_t B &= 0, \quad \dot{r}_x = \dot{r}_y = \dot{r}_z = 0.
\end{align*}
$$

As observed in [12], this system can be solved exactly and one can prove the Gauss law is preserved along
this step.

Hence, for $U^n = (f^n, E^n, B^n, r_*^n)$ given, one can update the numerical solution using the notation

$$
U^{n+1} = \varphi_{\Delta t}^{[1]}(U^n). \tag{40}
$$

4.2.2 Second subsystem

Regarding the second subsystem, using the notation $u = (r_x g_x, r_y g_y, r_z g_z) \in \mathbb{R}^3$, we consider

$$
\begin{align*}
\partial_t f + u \cdot \nabla_v f &= 0, \\
\partial_t E &= 0, \\
\partial_t B &= - \nabla_x \times u, \\
\dot{r}_x &= \frac{1}{2} \int g_x (-J_x + \partial_y B_z - \partial_z B_y) dx, \\
\dot{r}_y &= \frac{1}{2} \int g_y (-J_y - \partial_x B_z + \partial_z B_x) dx, \\
\dot{r}_z &= \frac{1}{2} \int g_z (-J_z + \partial_x B_y - \partial_y B_x) dx
\end{align*}
$$
Thanks to what has been done before, we propose the following time discretization

$$ \partial_t f + u^{n+1/2} \cdot \nabla_v f = 0, \qquad (42a) $$

$$ \begin{gathered} E^{n+1} = E^n, \\ B^{n+1} = B^n - \Delta t \nabla_x \times u^{n+1/2}, \end{gathered} \qquad (42b) $$

$$ r_x^{n+1} = r_x^n + \frac{\Delta t}{2} \int g_x(-J_x^{n+1/2} + \partial_y B_z^{n+1/2} - \partial_z B_y^{n+1/2}) dx, \qquad (42c) $$

$$ r_y^{n+1} = r_x^n + \frac{\Delta t}{2} \int g_y(-J_y^{n+1/2} - \partial_x B_y^{n+1/2} + \partial_z B_z^{n+1/2}) dx, \qquad (42d) $$

$$ r_z^{n+1} = r_z^n + \frac{\Delta t}{2} \int g_z(-J_z^{n+1/2} + \partial_x B_z^{n+1/2} - \partial_y B_x^{n+1/2}) dx, \qquad (42e) $$

where we used $r_{\star}^{n+1/2} = \frac{1}{2}(r_{\star}^n + r_{\star}^{n+1})$, $\star = x, y, z$ and $u_{\star}^{n+1/2} = r_{\star}^{n+1/2}g_{\star}$, $\star = x, y, z$. Regarding (42a), a semi-Lagrangian method will be used

$$ f^{n+1}(x, v) = f^n(x, v_x - \Delta t r_x^{n+1/2} g_x, v_y - \Delta t r_y^{n+1/2} g_y, v_z - \Delta t r_z^{n+1/2} g_z), \qquad (43) $$

for which a one-dimensional splitting can be performed since it is a linear transport equation. Now, we will manipulate the equations (18b), (42c), (42d) and (42e) to write down the solutions can be updated with an explicit cost. Indeed, we can see $J_{\star}^{n+1/2}$ and $B_{\star}^{n+1/2}$ can be expressed as linear function of $r_{\star}^{n+1/2}$ so that we will derive a linear system for $r_{\star}^{n+1/2}$.

First, let derive the relation for the current $J_x$. To do so, we multiply (43) by $v_x$ and integrate on $v$ (the same strategy is performed for the components $J_y$ and $J_z$) to get

$$ J_x^{n+1} = J_x^n + \Delta t r_x^{n+1/2} g_x \rho^n, $$

so that $J_x^{n+1/2} = (J_x^{n+1} + J_x^n)/2$ is given by

$$ J_x^{n+1/2} = J_x^n + \frac{\Delta t}{2} r_x^{n+1/2} g_x \rho^n. $$

Regarding now the Faraday equation (42b), we observe that $B^{n+1/2}$ can be expressed in terms of $r^{n+1/2}$ so that we get for $r_x^{n+1}$

$$ \begin{align} r_x^{n+1} &= r_x^n + \frac{\Delta t}{2} \int g_x(-J_x^{n+1/2} + \partial_y B_z^{n+1/2} - \partial_z B_y^{n+1/2}) dx, \tag{44} \\ &= r_x^n + \frac{\Delta t}{2} \int g_x(-J_x^n - \frac{\Delta t}{2} r_x^{n+1/2} g_x \rho^n) dx \nonumber \\ &\quad - \frac{\Delta t}{2} \int (\partial_y g_x)(B_z^n - \frac{\Delta t}{2} r_y^{n+1/2} \partial_x g_y + \frac{\Delta t}{2} r_x^{n+1/2} \partial_y g_x) dx \nonumber \\ &\quad + \frac{\Delta t}{2} \int (\partial_z g_x)(B_y^n + \frac{\Delta t}{2} r_z^{n+1/2} \partial_x g_z - \frac{\Delta t}{2} r_x^{n+1/2} \partial_z g_x) dx \nonumber \\ &= r_x^n - \frac{\Delta t}{2} \int g_x J_x^n dx - \frac{\Delta t}{2} \int ((\partial_y g_x) B_z^n - (\partial_z g_x) B_y^n) dx \nonumber \\ &\quad - \frac{\Delta t^2}{8} r_x^{n+1} \int (g_x^2 \rho^n + (\partial_y g_x)^2 + (\partial_z g_x)^2) dx \nonumber \\ &\quad - \frac{\Delta t^2}{8} r_x^n \int (g_x^2 \rho^n + (\partial_y g_x)^2 + (\partial_z g_x)^2) dx \nonumber \\ &\quad + \frac{\Delta t^2}{8} r_y^{n+1} \int \partial_y g_x \partial_x g_y dx - \frac{\Delta t^2}{8} r_y^n \int \partial_y g_x \partial_x g_y dx \nonumber \\ &\quad + \frac{\Delta t^2}{8} r_z^{n+1} \int \partial_z g_x \partial_x g_z dx - \frac{\Delta t^2}{8} r_z^n \int \partial_z g_x \partial_x g_z dx. \nonumber \end{align} $$
In this case, the terms $r_x^{n+1}$, $r_y^{n+1}$, $r_z^{n+1}$ are coupled so that we have to solve a 3x3 system

$$ 
\begin{pmatrix} 1 + \frac{\Delta t^2}{8} a_{xx} & -\frac{\Delta t^2}{8} a_{xy} & -\frac{\Delta t^2}{8} a_{xz} \\ -\frac{\Delta t^2}{8} a_{xy} & 1 + \frac{\Delta t^2}{8} a_{yy} & -\frac{\Delta t^2}{8} a_{yz} \\ -\frac{\Delta t^2}{8} a_{xz} & -\frac{\Delta t^2}{8} a_{zy} & 1 + \frac{\Delta t^2}{8} a_{zz} \end{pmatrix} \begin{pmatrix} r_x^{n+1} \\ r_y^{n+1} \\ r_z^{n+1} \end{pmatrix} = \begin{pmatrix} \frac{\Delta t^2}{8} a_{xx} & \frac{\Delta t^2}{8} a_{xy} & \frac{\Delta t^2}{8} a_{xz} \\ \frac{\Delta t^2}{8} a_{xy} & -\frac{\Delta t^2}{8} a_{yy} & \frac{\Delta t^2}{8} a_{yz} \\ \frac{\Delta t^2}{8} a_{xz} & \frac{\Delta t^2}{8} a_{zy} & -\frac{\Delta t^2}{8} a_{zz} \end{pmatrix} \begin{pmatrix} r_x^n \\ r_y^n \\ r_z^n \end{pmatrix} \\ + \frac{\Delta t}{2} \begin{pmatrix} \int (-g_x^n J_x^n - (\partial_y g_x) B_z^n + (\partial_z g_x) B_y^n) dx \\ \int (-g_y^n J_y^n + (\partial_x g_y) B_z^n - (\partial_z g_y) B_x^n) dx \\ \int (-g_z^n J_z^n - (\partial_x g_z) B_y^n + (\partial_y g_z) B_x^n) dx \end{pmatrix}, \quad (45) $$

with

$$ 
\begin{align} 
a_{xx} &= \int (g_x^2 \rho + (\partial_y g_x)^2 + (\partial_z g_x)^2) dx, \notag \\ 
a_{yy} &= \int (g_y^2 \rho + (\partial_x g_y)^2 + (\partial_z g_y)^2) dx, \notag \\ 
a_{zz} &= \int (g_z^2 \rho + (\partial_x g_z)^2 + (\partial_y g_z)^2) dx, \notag \\ 
a_{xy} &= \int (\partial_y g_x)(\partial_x g_y) dx, \quad a_{xz} = \int (\partial_z g_y)(\partial_x g_z) dx, \quad a_{yz} = \int (\partial_y g_z)(\partial_z g_y) dx. \tag{46} 
\end{align} $$

The invertibility of the coefficient matrix can be found in Appendix A.

Hence, for $U^n = (f^n, E^n, B^n, r_\star^n)$ given, one can update $r_\star^{n+1}$ from the 3x3 linear system (45), then compute $f^{n+1}$ from (43) and $B^{n+1}$ from (42b). These updates can be written as

$$ U^{n+1} = \varphi_{\Delta t}^{[2]}(U^n). \quad (47) $$

### 4.2.3 Third subsystem

Finally, the last subsystem can be written as

$$ 
\begin{align*} 
\partial_t f + (v \times B) \cdot \nabla_v f &= 0, \\ 
\partial_t E &= - \nabla_x \times B, \\ 
\partial_t B &= 0, \quad \dot{r}_x = \dot{r}_y = \dot{r}_z = 0. 
\end{align*} $$

This system can be solved exactly in time and the solution can be written as

$$ 
\begin{align*} 
f^{n+1}(x, v) &= f^n(x, e^{\Delta t \tilde{B}^n} v), \\ 
E^{n+1} &= E^n - \Delta t \nabla_x \times B^n \\ 
B^{n+1} &= B^n, \quad r_\star^{n+1} = r_\star^n, \star = x, y, z, 
\end{align*} $$

where $\tilde{B}$ is the skew symmetric matrix such that $\tilde{B}v = v \times B$. The update of $f^{n+1}$ requires a three-dimensional interpolation which may be costly from a computational point of view. However, as observed in [2, 1], see also [20], this advection can be split exactly as one-dimensional shears.

Hence, for $U^n = (f^n, E^n, B^n, r_\star^n)$ given, one can update the numerical solution using the notation

$$ U^{n+1} = \varphi_{\Delta t}^{[3]}(U^n). \quad (49) $$

## 4.3 Properties of the new splitting

Let consider the splitting presented in the previous subsections, that we can reformulate using the notations $\varphi_{\Delta t}^{[i]}$, $i = 1, 2, 3$ and $U^n = (f^n, E^n, B^n, r_\star^n)$

$$ U^{n+1} = \varphi_{\Delta t}^{[3]} \circ \varphi_{\Delta t}^{[2]} \circ \varphi_{\Delta t}^{[1]}(U^n), \quad (50) $$

where $\varphi_{\Delta t}^{[i]}$, $i = 1, 2, 3$ are given by (40), (47) and (49). One can prove the following result

**Proposition 6.** The scheme (50) preserves
• the Gauss law $\nabla_x \cdot E^n = \rho^n - m(0)$ for all $n \in \mathbb{N}$

• the following total energy $\tilde{\mathcal{E}}^n = \int \int \frac{|v|^2}{2} f^n(x, v) dxdv + (r_x^n)^2 + (r_y^n)^2 + (r_z^n)^2 + \frac{1}{2} \int |B^n|^2 dx$.

*Proof.* Here we use the notation $U^{n+1} = \varphi_{\Delta t}^{[3]} \circ \varphi_{\Delta t}^{[2]} \circ \varphi_{\Delta t}^{[1]}(U^n) = \varphi_{\Delta t}^{[3]} \circ \varphi_{\Delta t}^{[2]}(U^{(1)}) = \varphi_{\Delta t}^{[3]}(U^{(2)})$. Then, we observe that the kinetic energy is preserved along the first and third step so that $\mathcal{E}_{kin}^{(1)} = \mathcal{E}_{kin}^n$ and $\mathcal{E}_{kin}^{n+1} = \mathcal{E}_{kin}^{(2)}$. Moreover, $B, r_*$ are unchanged along these first and third steps, so we focus on the second step and as in the proof of Prop. 4, we consider the update from iteration $n$ to $n+1$ (instead of (2)) for this second step.

Let compute the kinetic energy $\mathcal{E}_{kin}^{n+1}$ after the second step

$$ 
\begin{aligned}
\mathcal{E}_{kin}^{n+1} &= \frac{1}{2} \iint |v|^2 f^{n+1}(x, v) dxdv = \frac{1}{2} \iint |v|^2 f^n(x, v - \Delta t u^{n+1/2}) dxdv \\
&= \frac{1}{2} \iint |w + \Delta t u^{n+1/2}|^2 f^n(x, w) dxdw \\
&= \mathcal{E}_{kin}^n + \Delta t \iint u^{n+1/2} \cdot \left(w + \frac{\Delta t}{2} u^{n+1/2}\right) f^n(x, w) dxdw \\
&= \mathcal{E}_{kin}^n + \Delta t \int u^{n+1/2} \cdot \left(J^n \frac{\Delta t}{2} \rho^n u^{n+1/2}\right) dx \\
&= \mathcal{E}_{kin}^n + \Delta t \int u^{n+1/2} \cdot J^{n+1/2} dx.
\end{aligned} 
\quad (51) $$

Let now consider the $r$ unknown and denote $r_i, i = x, y, z$ the components. Multiplying the scheme (42c), (42d) and (42e) on $r_i, i = x, y, z$ by $(r_i^{n+1} + r_i^n)$ leads to

$$ (r_i^{n+1})^2 - (r_i^n)^2 = -\Delta t \int r_i^{n+1/2} \int g_i J_i^{n+1/2} dx + \Delta t \int r_i^{n+1/2} \int g_i (\nabla \times B^{n+1/2})_i dx $$

and summing on $i = x, y, z$ leads to

$$ 
\begin{aligned}
\sum_{i=x,y,z} (r_i^{n+1})^2 - (r_i^n)^2 &= -\Delta t \int u^{n+1/2} \cdot J^{n+1/2} dx + \Delta t \int u^{n+1/2} \cdot (\nabla \times B^{n+1/2}) dx \\
&= -\Delta t \int u^{n+1/2} \cdot J^{n+1/2} dx + \Delta t \int (\nabla \times u^{n+1/2}) \cdot B^{n+1/2} dx.
\end{aligned} 
\quad (52) $$

Then, we multiply the equation (42b) on $B$ by $B^{n+1/2}$ and integrate on $x$ to get

$$ \frac{1}{2} \int (|B^{n+1}|^2 - |B^n|^2) dx = -\Delta t \int B^{n+1/2} \cdot (\nabla \times u^{n+1/2}) dx. \quad (53) $$

Finally, adding (51), (52) and (53) gives the total energy conservation for the second step and thus for the whole scheme. $\square$

# 5 Phase space approximation

In this section, some details are provided on the phase space approximation. A similar strategy as the one used in [12], but we would like to emphasize that the proposed splittings can be combined with other strategies.

## Space approximation
Here we used a Fourier spectral approximation in space which is a natural choice when periodic boundary conditions are considered. The computations are greatly simplified since spatial derivatives reduce to multiplication in Fourier space, in particular for the Maxwell equations.

Now, we give some details on the approximation of $\varphi_{\Delta t}^{[1]}$, the other steps being simple. Assuming $U^n = (f^n, E^n, B^n, r_\star^n)$ known, we will detail how to compute $U^{n+1} = \varphi_{\Delta t}^{[1]}(U^n)$. Let us recall the
equations to solve for this step

$$
\begin{align*}
\partial_t f + v \cdot \nabla_x f &= 0, \\
\partial_t E &= -J, \\
\partial_t B &= 0, \quad \dot{r}_x = \dot{r}_y = \dot{r}_z = 0.
\end{align*}
$$

We first observe $B^n, r_*^n$ are unchanged and let denote $\hat{f}_k^n, \hat{E}_k^n$ the Fourier coefficients of $f^n, E^n$ in space, $k \in \mathbb{Z}^3$. By using the Fourier transform of the first two equations, one obtains

$$
\begin{align*}
\frac{d\hat{f}_k}{dt} + v \cdot (ik)\hat{f}_k &= 0 \\
\frac{d\hat{E}_k}{dt} &= -\int v \hat{f}_k dv
\end{align*}
$$

The first equation can be solved exactly in time

$$
\hat{f}_k(t) = e^{-ik \cdot v(t-t_n)} \hat{f}_k^n,
$$

which can be inserted in the second equation integrated on $t \in [t_n, t_{n+1}]$ to get (for $k \neq (0,0,0)$)

$$
\begin{align*}
\hat{E}_k^{n+1} &= \hat{E}_k^n - \int_{t_n}^{t_{n+1}} \int v \hat{f}_k(s, v) dv ds \\
&= \hat{E}_k^n - \int \int_{t_n}^{t_{n+1}} v e^{-ik \cdot v(s-t_n)} \hat{f}_k^n(v) dv ds \\
&= \hat{E}_k^n + \int \frac{v}{ik \cdot v} (e^{-ik \cdot v \Delta t} - 1) \hat{f}_k^n(v) dv.
\end{align*}
$$

From this expression for $\hat{E}_k^{n+1}$, we easily deduce the Gauss law is propagated. Indeed, taking the scalar product of the previous relation with $ik$ leads to

$$
\begin{align*}
ik \cdot \hat{E}_k^{n+1} &= ik \cdot \hat{E}_k^n + \int \frac{ik \cdot v}{ik \cdot v} (e^{-ik \cdot v \Delta t} - 1) \hat{f}_k^n dv \\
&= \hat{\rho}_k^n - m(0) + \int \hat{f}_k^{n+1} dv - \int \hat{f}_k^n dv, \\
&= \int \hat{f}_k^{n+1} dv - m(0),
\end{align*}
$$

which means the discrete Poisson equation is satisfied.

Velocity approximation

Regarding the velocity direction, we essentially need to approximate integrals and one-dimensional lin-
ear transport equations (in the $\varphi_{\Delta t}^{[2]}$ and $\varphi_{\Delta t}^{[3]}$ steps). Regarding the integrals, we use standard trapezoidal
quadrature. Regarding the one-dimensional linear transport equations, a semi-Lagrangian method is used
with a high order Lagrange interpolation (of odd degree 2d + 1 for stability reasons [16, 8]). In more
details, to compute (43), a directional splitting (which is exact in this case) is used so that we have to
compute $f((v_x)_i - \Delta \text{tr}_x^{n+1/2} g_x, v_y, v_z)$ using high order interpolation, where $(v_x)_i$ denotes the grid point in
the $v_x$ direction. In practice, a Lagrange interpolation is utilized for the $v_x$ direction. Removing the $v_y, v_z$
direction to ease the reading, we introduce $j \in \mathbb{Z}$ such that $(v_x)_{i+j} \le (v_x)_i - \Delta \text{tr}_x^{n+1/2} g_x < (v_x)_{i+j+1}$
and $d \in \mathbb{Z}$ which defines the following centered $(2d+2)$-points stencil

$$
\{(v_x)_{i+j-d}, \dots, (v_x)_{i+j}, (v_x)_{i+j+1}, \dots, (v_x)_{i+j+1+d}\}.
$$

From the corresponding value of $f$ on this stencil, one can reconstruct a Lagrange polynomial $\mathcal{L}: v_x \in [(v_x)_{i+j}, (v_x)_{i+j+1}] \mapsto \mathcal{L}(v_x)$ of degree $(2d+1)$ which is used to approximate $f(v_x - \Delta \text{tr}_x^{n+1/2} g_x)$

$$
f((v_x)_i - \Delta \mathrm{tr}_x^{n+1/2} g_x) \approx \mathcal{L}((v_x)_i - \Delta \mathrm{tr}_x^{n+1/2} g_x).
$$

We refer to [15, 8, 7, 16] for more details.
# 6 Numerical results

In this part, we present some numerical results obtained by our approach we call. First, some tests on the Vlasov-Ampère model ($x, v \in \mathbb{R}$) and Vlasov-Maxwell VM1d-2v system ($x, v \in \mathbb{R} \times \mathbb{R}^2$) are shown. Let us remark the corresponding implementation has been performed using double-precision MATLAB or Julia programming. Second, we will discuss a JAX implementation [4] of the method using the VM1d-2v model, that uses single-precision.

## 6.1 Vlasov-Ampère

We first consider the Vlasov-Ampère model on the so-called Two Stream Instability (TSI). Hence, we consider the following initial condition

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

We compare three numerical schemes :

* VASAV which is the new scheme (9) described in Subsection 2.2.

* VA which is the scheme proposed in [19]

* VP which is the standard Hamiltonian splitting for the Vlasov-Poisson equations.

For all the schemes, we use a Fourier spectral method in space whereas a semi-Lagrangian scheme is used for the velocity transport term (the scheme SL-(2d + 1) from [15] with d = 9 is used). To illustrate the advantages of the new scheme, we consider the error of the Poisson equation $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$ for each iteration together with the relative total energy $\tilde{\mathcal{E}}(t) = \int \int \frac{v^2}{2} f(t, x, v) dxdv + r^2(t)$. For VA and VP, we also plot the quantity $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$ whereas the total energy $\mathcal{E}(t) = \int \int \frac{v^2}{2} f(t, x, v) dxdv + \frac{1}{2} \int E^2(t, x) dx$ is considered. Let us recall the VA scheme preserves the total energy $\mathcal{E}(t)$ but does not preserve the charge. For VP, the scheme preserved the Poisson equation by construction since the Poisson equation is resolved, but the total energy $\mathcal{E}(t)$ is not ensured.

The numerical parameters are chosen as follows $N_x = N_v = 128$ and $\Delta t = 0.5$. In Figure 1 are displayed the time evolution of $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$ and of $\tilde{\mathcal{E}}(t)$ for VASAV and $\mathcal{E}(t)$ for VA and VP. One can observe the new method preserves both the Gauss law and $\tilde{\mathcal{E}}(t)$ (up to machine accuracy) whereas VA does not preserve the Gauss law and VP does not preserve the total energy $\mathcal{E}(t)$.

<figure><img src="image_4.png" /><figcaption>Figure 1: Left: time evolution of the error on the Poisson equation $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$. Comparison between method VA from [19] and VASAV. Right: time evolution of the error on the total energy. Comparison between method VP and VASAV.</figcaption></figure>

On Figure 2, the error between the scalar auxiliary variable $r^2$ and the electric energy $\int E^2 dx$ is plotted as a function of the time step (semi-log scale). The time step are chosen in the range $\Delta t \in \{10^{-p}, p = 1, \dots, 6\}$. As expected, the second order accuracy is observed. This illustrates the effect of the SAV formulation.
<figure><img src="image_5.png" /><figcaption>Figure 2: Left: SAV error between $r^2$ and $\int E^2(t, x)dx$ for different time steps $\Delta t$ (semi-log scale). Right: Phase-space solution at the final time $t = 100$ using VASAV.</figcaption></figure>

## 6.2 Weibel test

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

$$f(t = 0, x, v_x, v_y) = \frac{1}{2\pi\sigma_1\sigma_2} \exp\left(-\frac{1}{2}\left(\frac{v_x^2}{\sigma_1^2} + \frac{v_y^2}{\sigma_2^2}\right)\right), \quad x \in [0, 2\pi/k), (v_x, v_y) \in \mathbb{R}^2,$$ $$B(t = 0, x) = \beta \cos(kx), \quad E_y(t = 0, x) = 0,$$

and $E_x(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$, $\beta = -10^{-4}$ for our test. This choice of initial condition, induces a Weibel type instability with a growth rate of 0.02784 (obtained by solving the dispersion relation, see [26]).

For the numerical simulations up to a final time $T = 300$, the numerical parameters are chosen as follows: $N_x = 16$, $N_{v_x} = N_{v_y} = 47$, the velocity domain is $(v_x, v_y) \in [-0.2, 0.2] \times [-0.5, 0.5]$ and $\Delta t = 1/16$.

As for the Vlasov-Ampère case, we compare three numerical schemes:

* VMSAV which is the new scheme (50) described in Subsection 3.2

* ECSL which is the scheme proposed in [19]

* VM-split which is the Hamiltonian splitting for the Vlasov-Maxwell equations [12].

For all the schemes, we use a Fourier spectral method in space whereas a semi-Lagrangian scheme is used for the velocity transport term (again the SL-(2d + 1) is used). To illustrate the advantages of the new scheme, we consider the error of the Gauss law $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$ for each iteration together with the relative energy $\tilde{\mathcal{E}}^n$ for VMSAV and $\mathcal{E}^n$ for ECSL and VM-split.

In Figures 3, we plot the time evolution of the three electro-magnetic energies (semi-log scale), namely $\|B(t)\|_{L_x^2}$, $\|E_x(t)\|_{L_x^2}$ and $\|E_y(t)\|_{L_x^2}$. The three solvers gives very similar results, and the growth rate calculated for $\|B(t)\|_{L_x^2}$ is about 0.054 which is very close to the one obtained from the linear theory ($2 \times 0.02784$). The results are also in very good agreement with ones obtained in the literature see [17, 12, 14]. Then in Figure 4-left, we plot the time evolution of $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$ for the new method VMSAV and ECSL (this quantity is preserved by VM-split). One can observe the ECSL solver does not preserve exactly the Gauss law. In Figure 4-right, we plot the time history of $\tilde{\mathcal{E}}$ for VMSAV and of $\mathcal{E}$ for VM-split ($\mathcal{E}$ is preserved exactly for ECSL). As expected, the total energy $\tilde{\mathcal{E}}$ is preserved up to machine accuracy for VMSAV whereas VM-split does not preserve $\mathcal{E}$.
<figure><img src="image_6.png" /><figcaption>Figure 3: Time evolution of the magnetic and electric energies (semi-log scale) for the 3 VM solvers: VMSAV (left), ECSL (middle), VM-split (right).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Left: time evolution of the error on the Poisson equation $\|\partial_x E(t, \cdot) - (\int f(t, \cdot, v)dv - 1)\|_{L_x^2}$. Comparison between method ECSL from [19] and VMSAV. Right: time evolution of the error on the total energy. Comparison between the splitting method from [12] and VMSAV.</figcaption></figure>

## 6.3 Landau-Weibel test

We consider a second test which combines an instability for the $B$ and $E_y$ field and a damping for the $E_x$ field. This can be obtained using the following initial condition

$$f(t = 0, x, v_x, v_y) = \frac{1}{2\pi\sigma_1\sigma_2} \exp\left(-\frac{1}{2}\left(\frac{v_x^2}{\sigma_1^2} + \frac{v_y^2}{\sigma_2^2}\right)\right) (1 + \alpha \cos(kx), \quad x \in [0, 2\pi/k),$$

$$B(t = 0, x) = \beta \cos(kx), \quad E_y(t = 0, x) = 0,$$

and $E_x(x, t = 0)$ is imposed from the Poisson equation. We choose the parameters $\sigma_1 = 1$, $\sigma_2 = \sqrt{2}\sigma_1$, $k = 0.5$, $\alpha = 10^{-3}$ and $\beta = -10^{-6}$ for our test.

The numerical parameters are $N_x = N_{v_x} = N_{v_y} = 47$ and $\Delta t = 0.125$ for a final time $T = 50$. From the linear theory, one can check the longitudinal electric field $E_x$ enjoys a Landau damping with a rate $\gamma_L = -0.153 < 0$ whereas the transverse fields $B$, $E_y$ are unstable. These two phenomena can be observed in Figure 5 where the time evolution of $\|B\|_{L_x^2}^2$, $\|E_y\|_{L_x^2}^2$ and $\|E_x\|_{L_x^2}^2$ are displayed for the three solvers VMSAV, ECSL and VM-split.

## 6.4 JAX Acceleration

We reimplemented the VMSAV (one dimension in space and two dimensions in velocity) solver in JAX [4] and leveraged XLA compilation [9] to target NVIDIA Quadro RTX 6000 GPUs. This approach yields a substantial performance boost by taking advantage of JAX’s just-in-time compilation and GPU-optimized kernels. On the hand, due to memory constraints, we use JAX’s default single-precision arithmetic rather than double precision.
<figure><img src="image_8.png" /><figcaption>Figure 5: Time evolution of the magnetic and electric energies (semi-log scale) for the 3 VM solvers: VMSAV (left), ECSL (middle), VM-split (right).</figcaption></figure>

Before discussing the performances, let us precise the main differences of the JAX implementation compared to the scheme presented in Section 3:

* the first step is the same and Fourier approximation is used in space.

* for the semi-Lagrangian method in the second step, the one-dimensional interpolations are performed with JAX's built-in linear interpolator `jax numpy interp`;

* for the third step, the two-dimensional cubic B-spline interpolations in velocity space are carried out by our custom, JIT-compiled function.

To evaluate performance, we executed both the Weibel instability test (Section 6.2, final time $T = 300$) and the Landau-Weibel test (Section 6.3, final time $T = 50$) 100 times across a range of time-step sizes $\Delta t$ and grid resolutions (denoted as $N_x \times N_{v_x} \times N_{v_y}$), and recorded the average runtime. The results are presented in Tables 1 and 2.

<table><caption>Table 1: Average runtime for the Weibel instability over 100 runs.</caption><thead><tr><th>$N_x \times N_{v_x} \times N_{v_y}$</th><th>8 × 23 × 23</th><th>16 × 47 × 47</th><th>32 × 95 × 95</th></tr></thead><tbody><tr><th scope="row">0.12500</th><td>1.3809s</td><td>1.9298s</td><td>6.4328s</td></tr><tr><th scope="row">0.06250</th><td>2.7623s</td><td>3.8784s</td><td>12.8695s</td></tr><tr><th scope="row">0.03125</th><td>5.5249s</td><td>7.7480s</td><td>25.7391s</td></tr><tr><th scope="row">0.01562</th><td>11.0874s</td><td>15.4943s</td><td>51.3865s</td></tr></tbody></table>

<table><caption>Table 2: Average runtime for the Landau-Weibel test over 100 runs.</caption><thead><tr><th>$N_x \times N_{v_x} \times N_{v_y}$</th><th>23 × 23 × 23</th><th>47 × 47 × 47</th><th>95 × 95 × 95</th><th>191 × 191 × 191</th></tr></thead><tbody><tr><th scope="row">0.25000</th><td>0.1111s</td><td>0.2655s</td><td>1.8103s</td><td>15.9224s</td></tr><tr><th scope="row">0.12500</th><td>0.2219s</td><td>0.5315s</td><td>3.6211s</td><td>31.8685s</td></tr><tr><th scope="row">0.06250</th><td>0.4426s</td><td>1.0643s</td><td>7.2388s</td><td>63.7605s</td></tr><tr><th scope="row">0.03125</th><td>0.8836s</td><td>2.1203s</td><td>14.4742s</td><td>127.6236s</td></tr></tbody></table>

Both tables indicate that runtime scales approximately linearly with $1/\Delta t$. Table 1 shows that doubling the spatial resolution increases runtime by a factor of roughly 1–2 up to 3–4. Similarly, Table 2 confirms near-linear scaling with time resolution, while scaling in the spatial and velocity domains deviates from linearity. It is likely due to the semi-Lagrangian scheme and interpolation step it requires. Overall, these results demonstrate that our JAX-XLA implementation delivers excellent GPU-accelerated performance with predictable scaling behavior in time.

Finally, we present the results obtained with our JAX-XLA implementation for the Weibel test case (see Section 6.2). Using a refined mesh of $N_x \times N_{v_x} \times N_{v_y} = 64 \times 191 \times 191$ and a time step $\Delta t = 0.01562$, Figure 6 (left) shows the time evolution of the electromagnetic energies on a semi-logarithmic scale, highlighting the instabilities in $\|E_x\|_{L^2}$, $\|E_y\|_{L^2}$, and $\|B\|_{L^2}$. The right panel of Figure 6 shows the
<figure><img src="image_9.png" /><figcaption>Figure 6: Left: time evolution of electromagnetic energies $||E_x||_{L^2}$, $||E_y||_{L^2}$, $||B||_{L^2}$ in semi-log scale. Right: time evolution of electromagnetic energies $||E_x||_{L^2}$, $||E_y||_{L^2}$, $||B||_{L^2}$ and total energy $\tilde{\mathcal{E}}$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 7: Left: Time evolution of $\tilde{\mathcal{E}}(t) - \tilde{\mathcal{E}}(0)$. Right: 3D scatter of $(v_x, v_y) \mapsto f(t = 300, x, v_x, v_y)$ at four x-locations $(x = 0, 1.65, 3.30, 4.95)$.</figcaption></figure>

electromagnetic energies alongside the kinetic energy (denoted by $K$), the total modified energy ($\tilde{\mathcal{E}}$) and the change in total modified energy, $\Delta\tilde{\mathcal{E}} = \tilde{\mathcal{E}}(t) - \tilde{\mathcal{E}}(0)$, on a linear scale. This plot showcases the energy-exchange dynamics: part of the kinetic energy is converted into magnetic energy, whereas the variation in $\tilde{\mathcal{E}}$ remains vanishingly small. Figure 7 (left) presents the relative change in the total energy $\tilde{\mathcal{E}}$, demonstrating preservation close to single-precision accuracy. The right panel shows the distribution function $f(t = 300, x, v_x, v_y)$ at spatial locations $x = 0, 1.65, 3.30, 4.95$, revealing clear deviations from the initial equilibrium.

# 7 Conclusion

In this work, a novel scheme is proposed for the approximation of the Vlasov-Ampère (VA) and Vlasov-Maxwell (VM) equations. The numerical scheme is based on a reformulation of the VA or VM model where an additional scalar variable is introduced, inspired by the SAV approach. This scheme enjoys the following properties: (i) it preserves the Gauss law without solving it; (ii) there is no stringent CFL condition thanks to the use of semi-Lagrangian scheme; (iv) the method is explicit so its cost is the same as the one of an explicit solver. Some tests are conducted on the VA system (1d-1v) and on the VM1d-2v model which illustrate the good behavior of the scheme compared to the solvers of the literature that only preserve the total energy or the Gauss law.

Some extensions of this work would be the following: different space approximations will be considered (to avoid the use of Fourier approximation which can not be extended to complex geometries/meshes), alternative splitting inspired from [20] may be useful to overcome the CFL condition coming from the
Maxwell equations, or the extension of the VMSAV formulation for PIC methods.

A Invertibility of a matrix

In this Appendix, we study the invertibility of the matrix involved in (45).

**Lemme 1.** Let $\rho(x) \ge 0$ be a nonnegative weight function and let $g_x, g_y, g_z$ be functions such that the following integrals are finite:

$$
\begin{align*}
a_{xx} &= \int (g_x^2 \rho + (\partial_y g_x)^2 + (\partial_z g_x)^2) dx, & a_{xy} &= \int (\partial_y g_x)(\partial_x g_y) dx, \\
a_{yy} &= \int (g_y^2 \rho + (\partial_x g_y)^2 + (\partial_z g_y)^2) dx, & a_{xz} &= \int (\partial_z g_x)(\partial_x g_z) dx, \\
a_{zz} &= \int (g_z^2 \rho + (\partial_x g_z)^2 + (\partial_y g_z)^2) dx, & a_{yz} &= \int (\partial_z g_y)(\partial_y g_z) dx.
\end{align*}
$$

Define the matrix

$$
M = \begin{pmatrix}
1 + \frac{\Delta t^2}{8} a_{xx} & -\frac{\Delta t^2}{8} a_{xy} & -\frac{\Delta t^2}{8} a_{xz} \\
-\frac{\Delta t^2}{8} a_{xy} & 1 + \frac{\Delta t^2}{8} a_{yy} & -\frac{\Delta t^2}{8} a_{yz} \\
-\frac{\Delta t^2}{8} a_{xz} & -\frac{\Delta t^2}{8} a_{yz} & 1 + \frac{\Delta t^2}{8} a_{zz}
\end{pmatrix}.
$$

Then for any $\Delta t > 0$, $M$ is symmetric positive definite and hence invertible.

Proof. Let $c = (c_x, c_y, c_z)^T \in \mathbb{R}^3$ be arbitrary and nonzero, and set $\alpha = \Delta t^2/8$. Since $M$ is real symmetric, it suffices to show

$$
c^T M c = c^T I c + \alpha c^T B c = \|c\|^2 + \alpha c^T B c > 0,
$$

where the auxiliary matrix $B$ has entries $B_{ii} = a_{ii}$ and $B_{ij} = -a_{ij}$ for $i \neq j$.

First we compute $c^T B c$. By definition,

$$
\begin{align*}
c^T B c &= c_x^2 a_{xx} + c_y^2 a_{yy} + c_z^2 a_{zz} - 2c_x c_y a_{xy} - 2c_x c_z a_{xz} - 2c_y c_z a_{yz} \\
&= \int \left[ c_x^2 (g_x^2 \rho + (\partial_y g_x)^2 + (\partial_z g_x)^2) + c_y^2 (g_y^2 \rho + (\partial_x g_y)^2 + (\partial_z g_y)^2) \right. \\
&\qquad \left. + c_z^2 (g_z^2 \rho + (\partial_x g_z)^2 + (\partial_y g_z)^2) - 2c_x c_y (\partial_y g_x)(\partial_x g_y) - 2c_x c_z (\partial_z g_x)(\partial_x g_z) - 2c_y c_z (\partial_z g_y)(\partial_y g_z) \right] dx.
\end{align*}
$$

Next we regroup the integrand into sums of squares:

$$
\begin{align*}
c_x^2 g_x^2 \rho + c_y^2 g_y^2 \rho + c_z^2 g_z^2 \rho &= (c_x g_x \sqrt{\rho})^2 + (c_y g_y \sqrt{\rho})^2 + (c_z g_z \sqrt{\rho})^2, \\
c_x^2 (\partial_y g_x)^2 + c_y^2 (\partial_x g_y)^2 - 2c_x c_y (\partial_y g_x) (\partial_x g_y) &= [c_x \partial_y g_x - c_y \partial_x g_y]^2, \\
c_x^2 (\partial_z g_x)^2 + c_z^2 (\partial_x g_z)^2 - 2c_x c_z (\partial_x g_z) (\partial_x g_z) &= [c_x \partial_z g_x - c_z \partial_x g_z]^2, \\
c_y^2 (\partial_z g_y)^2 + c_z^2 (\partial_x g_z)^2 - 2c_y c_z (\partial_x g_z) (\partial_x g_z) &= [c_y \partial_z g_y - c_z \partial_x g_z]^2.
\end{align*}
$$

Therefore

$$
c^T B c = \int \left[ (c_x g_x \sqrt{\rho})^2 + (c_y g_y \sqrt{\rho})^2 + (c_z g_z \sqrt{\rho})^2 + [c_x \partial_y g_x - c_y \partial_x g_y]^2 \right. \\
\left. + [c_x \partial_z g_x - c_z \partial_x g_z]^2 + [c_y \partial_z g_y - c_z \partial_y g_z]^2 \right] dx \ge 0.
$$

Since $\|c\|^2 > 0$ for $c \neq 0$, we conclude

$$
c^T M c = \|c\|^2 + \alpha c^T B c \geq \|c\|^2 > 0.
$$

Hence $M$ is positive definite and invertible for any $\Delta t > 0$. $\square$

<figure><img src="image_11.png" /></figure>
**Acknowledgement.** This work has been carried out within the framework of the EUROfusion Con-
sortium, funded by the European Union via the Euratom Research and Training Programme (grant
agreement No 101052200 EUROfusion). Views and opinions expressed are however those of the authors
only and do not necessarily reflect those of the European Union or the European Commission. NC has
been supported by the French Federation for Magnetic Fusion Studies (FR-FCM).

References

[1] J. Bernier, F. Casas, and N. Crouseilles. Splitting methods for rotation: applications to Vlasov equations. *SIAM J. Sc. Comput.*, 42:666–697, 2020.

[2] J. Bernier, N. Crouseilles, and Y. Li. Exact splitting methods for kinetic and Schrödinger equations. *J. Sci. Comput.*, 86:1–35, 2021.

[3] S. Bilbao, M. Ducceschi, and F. Zama. Explicit exactly energy-conserving methods for Hamiltonian systems. *J. Comput. Phys.*, 472:111697, 2023.

[4] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang. JAX: composable transformations of Python+NumPy programs, 2018.

[5] M. Campos Pinto, S. Jund, S. Salmon, and E. Sonnendrücker. Charge-conserving FEM-PIC schemes on general grids. *C. R. Méc.*, 342(10):570–582, 2014.

[6] M. Campos Pinto and V. Pagès. A semi-implicit electromagnetic FEM-PIC scheme with exact energy and charge conservation. *J. Comput. Phys.*, 453:110912, 2022.

[7] F. Casas, N. Crouseilles, E. Faou, and M. Mehrenberger. High-order Hamiltonian splitting for Vlasov-Poisson equations. *Numer. Math.*, 135:769–801, 2017.

[8] F. Charles, B. Després, and M. Mehrenberger. Enhanced convergence estimates for semi-Lagrangian schemes. Application to the Vlasov-Poisson equation. *SIAM J. on Numer. Anal.*, 51:840–863, 2013.

[9] Y. Chen, Z. Wan, Y. Zhuang, N. Liu, D. Lo, and X. Yang. Understanding the oss communities of deep learning frameworks: A comparative case study of pytorch and tensorflow. *ACM Trans. Softw. Eng. Methodol.*, 34(3), 2025.

[10] C. Cheng and G. Knorr. The integration of the Vlasov equation in configuration space. *J. Comput. Phys.*, 22:330–351, 1976.

[11] Y. Cheng, A. J. Christlieb, and X. Zhong. Energy-conserving discontinuous Galerkin methods for the Vlasov-Maxwell system. *J. Comput. Phys.*, 279:145–173, 2014.

[12] N. Crouseilles, L. Einkemmer, and E. Faou. Hamiltonian splitting for the Vlasov-Maxwell equations. *J. Comput. Phys.*, 283:224–240, 2015.

[13] N. Crouseilles, L. Einkemmer, and E. Faou. An asymptotic preserving scheme for the relativistic Vlasov-Maxwell equations in the classical limit. *Comput. Phys. Commun.*, 209:13–26, 2016.

[14] N. Crouseilles and X. Hong. Exponential DG methods for Vlasov equations. *J. Comput. Phys.*, 498:112682, 2024.

[15] N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker. Conservative semi-Lagrangian schemes for Vlasov equations. *J. Comput. Phys.*, 229:1927–1953, 2010.

[16] B. Després. Uniform asymptotic stability of Strang’s explicit compact schemes for linear advection. *SIAM J. on Numer. Anal.*, 47(5):3956–3976, 2009.

[17] Y. Li, Y. He, Y. Sun, J. Niesen, H. Qin, and J. Liu. Solving the Vlasov-Maxwell equations using Hamiltonian splitting. *J. Comput. Phys.*, 396:381–399, 2019.
[18] H. Liu, X. Cai, Y. Cao, and G. Lapenta. Conservative semi-Lagrangian kinetic scheme coupled with implicit finite element field solver for multidimensional Vlasov Maxwell system. *Commun. Nonlinear Sci. Numer. Simul.*, 102:105941, 2021.

[19] H. Liu, X. Cai, Y. Cao, and G. Lapenta. An efficient energy conserving semi-Lagrangian kinetic scheme for the Vlasov-Ampère system. *J. Comput. Phys.*, 492:112412, 2023.

[20] H. Liu, C. Lu, G. Xia, R. Keppens, and G. Lapenta. An efficient energy conserving semi-Lagrangian kinetic scheme for the Vlasov-Maxwell system. *J. Comput. Phys.*, 2025.

[21] J. Qiu and A. Christlieb. A conservative high order semi-Lagrangian WENO method for the Vlasov equation. *J. Comput. Phys.*, 229:1130–1149, 2010.

[22] J. A. Rossmanith and D. C. Seal. A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations. *J. Comput. Phys.*, 230:6203–6232, 2011.

[23] J. Shen, J. Xu, and J. Yang. The scalar auxiliary variable (SAV) approach for gradient flows. *J. Comput. Phys.*, 353:407–416, 2018.

[24] N. J. Sircombe and T. D. Arber. Valis: A split-conservative scheme for the relativistic 2d Vlasov-Maxwell system. *J. Comput. Phys.*, 228:4773–4788, 2009.

[25] A. Stock, J. Neudorfer, R. Schneider, C. Altmann, and C.-D. Munz. Investigation of the purely hyperbolic Maxwell system for divergence cleaning in discontinuous Galerkin based particle-in-cell methods. *Coupled Problems 2011 IV International Conference on Computational Methods for Coupled Problems in Science and Engineering*, 2011.

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