# Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field

Francis Filbet, Luis Miguel Miguel Rodrigues

► To cite this version:

Francis Filbet, Luis Miguel Miguel Rodrigues. Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field. SIAM Journal on Numerical Analysis, 2016, 54 (2), pp.1120-1146. 10.1137/15M104952X. hal-01232404

HAL Id: hal-01232404

https://hal.science/hal-01232404

Submitted on 23 Nov 2015

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# ASYMPTOTICALLY STABLE PARTICLE-IN-CELL METHODS FOR THE VLASOV-POISSON SYSTEM WITH A STRONG EXTERNAL MAGNETIC FIELD

FRANCIS FILBET, LUIS MIGUEL RODRIGUES

**ABSTRACT.** This paper deals with the numerical resolution of the Vlasov-Poisson system with a strong external magnetic field by Particle-In-Cell (PIC) methods. In this regime, classical PIC methods are subject to stability constraints on the time and space steps related to the small Larmor radius and plasma frequency. Here, we propose an asymptotic-preserving PIC scheme which is not subjected to these limitations. Our approach is based on first and higher order semi-implicit numerical schemes already validated on dissipative systems [6]. Additionally, when the magnitude of the external magnetic field becomes large, this method provides a consistent PIC discretization of the guiding-center equation, that is, incompressible Euler equation in vorticity form. We propose several numerical experiments which provide a solid validation of the method and its underlying concepts.

**KEYWORDS.** High order time discretization; Vlasov-Poisson system; Guiding-centre model; Particle methods.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. A brief review of particle methods</td><td>4</td></tr><tr><td>3. A particle method for Vlasov-Poisson system with a strong magnetic field</td><td>6</td></tr><tr><td>3.1. A first order semi-implicit scheme</td><td>7</td></tr><tr><td>3.2. Second order semi-implicit Runge-Kutta schemes</td><td>8</td></tr><tr><td>3.3. Third order semi-implicit Runge-Kutta schemes</td><td>10</td></tr><tr><td>4. Analysis of the first-order semi-implicit scheme</td><td>12</td></tr><tr><td>5. Numerical simulations</td><td>15</td></tr><tr><td>5.1. One single particle motion</td><td>15</td></tr><tr><td>5.2. Diocotron instability</td><td>18</td></tr><tr><td>6. Conclusion and perspective</td><td>21</td></tr><tr><td>Acknowledgements</td><td>23</td></tr><tr><td>References</td><td>23</td></tr></table>

## 1. INTRODUCTION

Magnetized plasmas are encountered in a wide variety of astrophysical situations, but also in magnetic fusion devices such as tokamaks, where a large external magnetic field needs to be applied in order to keep the particles on the desired tracks.

In Particle-In-Cell (PIC) simulations of such devices, this large external magnetic field obviously needs to be taken into account when pushing the particles. However, due to the magnitude of the concerned field this often adds a new time scale to the simulation and thus a stringent restriction on the time step. In order to get rid of this additional time scale, we would like to find approximate equations, where only the gross behavior implied by the external field would
be retained and which could be used in a numerical simulation. In the simplest situation, the trajectory of a particle in a constant magnetic field **B** is a helicoid along the magnetic field lines with a radius proportional to the inverse of the magnitude of **B**. Hence, when this field becomes very large the particle gets trapped along the magnetic field lines. However due to the fast oscillations around the apparent trajectory, its apparent velocity is smaller than the actual one. This result has been known for some time as the guiding center approximation, and the link between the real and the apparent velocity is well known in terms of **B**.

Here, we consider a plasma constituted of a large number of charged particles, which is described by the Vlasov equation coupled with the Maxwell or Poisson equations to compute the self-consistent fields. It describes the evolution of a system of particles under the effects of external and self-consistent fields. The unknown $f(t, \mathbf{x}, \mathbf{v})$, depending on the time $t$, the position $\mathbf{x}$, and the velocity $\mathbf{v}$, represents the distribution of particles in phase space for each species with $(\mathbf{x}, \mathbf{v}) \in \mathbb{R}^d \times \mathbb{R}^d$, $d = 1, .., 3$. Its behaviour is given by the Vlasov equation,

$$ (1.1) \qquad \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f + \mathbf{F}(t, \mathbf{x}, \mathbf{v}) \cdot \nabla_{\mathbf{v}} f = 0, $$

where the force field $\mathbf{F}(t, \mathbf{x}, \mathbf{v})$ is coupled with the distribution function $f$ giving a nonlinear system. We first define $\rho(t, \mathbf{x})$ the charge density and $\mathbf{J}(t, \mathbf{x})$ the current density which are given by

$$ \rho(t, \mathbf{x}) = q \int_{\mathbb{R}^d} f(t, \mathbf{x}, \mathbf{v}) d\mathbf{v}, \quad \mathbf{J}(t, \mathbf{x}) = q \int_{\mathbb{R}^d} \mathbf{v} f(t, \mathbf{x}, \mathbf{v}) d\mathbf{v}, $$

where $q$ is the elementary charge. For the Vlasov-Poisson model

$$ (1.2) \qquad \mathbf{F}(t, \mathbf{x}, \mathbf{v}) = \frac{q}{m} \mathbf{E}(t, \mathbf{x}), \quad \mathbf{E}(t, \mathbf{x}) = -\nabla_{\mathbf{x}}\phi(t, \mathbf{x}), \quad -\Delta_{\mathbf{x}}\phi = \frac{\rho}{\varepsilon_0}, $$

where $m$ represents the mass of one particle. On the other hand for the Vlasov-Maxwell model, we have

$$ \mathbf{F}(t, \mathbf{x}, \mathbf{v}) = \frac{q}{m} (\mathbf{E}(t, \mathbf{x}) + \mathbf{v} \wedge \mathbf{B}(t, \mathbf{x})), $$

and $\mathbf{E}, \mathbf{B}$ are solutions of the Maxwell equations

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{E}}{\partial t} - c^2 \nabla \times \mathbf{B} = -\frac{\mathbf{J}}{\varepsilon_0}, \\[1.5ex] \displaystyle \frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{E} = 0, \\[1.5ex] \displaystyle \nabla \cdot \mathbf{E} = \frac{\rho}{\varepsilon_0}, \quad \nabla \cdot \mathbf{B} = 0, \end{array} \right. $$

with the compatibility condition

$$ \frac{\partial \rho}{\partial t} + \mathrm{div}_{\mathbf{x}} \mathbf{J} = 0, $$

which is verified by the solutions of the Vlasov equation.

Here we will consider an intermediate model where the magnetic field is given,

$$ \mathbf{B}(t, \mathbf{x}) = \frac{1}{\varepsilon} \mathbf{B}_{\mathrm{ext}}(t, \mathbf{x}_{\perp}), $$

with $\varepsilon > 0$ and we focus on the long time behavior of the plasma in the orthogonal plane to the external magnetic field, that is the two dimensional Vlasov-Poisson system with an external
strong magnetic field

$$ (1.3) \qquad \left\{ \begin{array}{l} \varepsilon \frac{\partial f}{\partial t} + \mathbf{v}_{\perp} \cdot \nabla_{\mathbf{x}_{\perp}} f + \left( \mathbf{E}(t, \mathbf{x}_{\perp}) + \frac{1}{\varepsilon} \mathbf{v}_{\perp} \wedge \mathbf{B}_{\text{ext}}(t, \mathbf{x}_{\perp}) \right) \cdot \nabla_{\mathbf{v}_{\perp}} f = 0, \quad (\mathbf{x}_{\perp}, \mathbf{v}_{\perp}) \in \mathbb{R}^4, \\ E = -\nabla_{\mathbf{x}_{\perp}} \phi, \quad -\Delta_{\mathbf{x}_{\perp}} \phi = \rho, \quad \mathbf{x}_{\perp} \in \mathbb{R}^2. \end{array} \right. $$

Here, for simplicity we set all physical constants to one and consider that $\varepsilon > 0$ is a small parameter related to the ratio between the reciprocal Larmor frequency and the advection time scale. The term $\varepsilon$ in front of the time derivative of $f$ stands for the fact that we want to approximate the solution for large time.

We want to construct numerical solutions to the Vlasov-Poisson system (1.3) by particle methods (see [5]), which consist in approximating the distribution function by a finite number of macro-particles. The trajectories of these particles are computed from the characteristic curves corresponding to the Vlasov equation

$$ (1.4) \qquad \left\{ \begin{array}{l} \varepsilon \frac{d\mathbf{X}}{dt} = \mathbf{V}, \\ \varepsilon \frac{d\mathbf{V}}{dt} = \frac{1}{\varepsilon} \mathbf{V} \wedge \mathbf{B}_{\text{ext}}(t, \mathbf{X}) + \mathbf{E}(t, \mathbf{X}), \\ \mathbf{X}(t^0) = \mathbf{x}_{\perp}^0, \quad \mathbf{V}(t^0) = \mathbf{v}_{\perp}^0, \end{array} \right. $$

where the electric field is computed from a discretization of the Poisson equation in (1.3) on a mesh of the physical space.

The main purpose of this work is the construction of efficient numerical methods for stiff transport equations of type (1.3) in the limit $\varepsilon \to 0$. Indeed, setting

$$ \mathbf{Z} = \frac{\mathbf{V}}{\varepsilon} - \frac{\mathbf{E} \wedge \mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2}, $$

the system (1.4) can be re-written for $(\mathbf{X}, \mathbf{Z})$ as

$$ (1.5) \qquad \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{X}}{dt} = \frac{\mathbf{E} \wedge \mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2} + \mathbf{Z}, \\ \displaystyle \frac{d\mathbf{Z}}{dt} = -\frac{1}{\varepsilon^2} \mathbf{B}_{\text{ext}} \wedge \mathbf{Z} - \frac{d}{dt} \left( \frac{\mathbf{E} \wedge \mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2} \right), \\ \mathbf{X}(t^0) = \mathbf{x}_{\perp}^0, \quad \mathbf{Z}(t^0) = \frac{1}{\varepsilon} \mathbf{v}_{\perp}^0 - \frac{\mathbf{E} \wedge \mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2} (t^0, \mathbf{x}_{\perp}^0). \end{array} \right. $$

Therefore, we denote by $(\mathbf{X}^{\varepsilon}, \mathbf{Z}^{\varepsilon})$ the solution to (1.5), and under some classical smoothness assumptions on the electromagnetic fields $(\mathbf{E}, \mathbf{B}_{\text{ext}})$, it is well-known at least when $\mathbf{v}_{\perp}^0 = 0$ or $\mathbf{B}_{\text{ext}}$ is homogeneous that $(\mathbf{Z}^{\varepsilon})_{\varepsilon>0}$ converges weakly to zero when $\varepsilon \to 0$, and $\mathbf{X}^{\varepsilon} \to \mathbf{Y}$, where $\mathbf{Y}$ corresponds to the guiding center approximation

$$ (1.6) \qquad \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{Y}}{dt} = \frac{\mathbf{E} \wedge \mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2} (t, \mathbf{Y}), \\ \mathbf{Y}(t^0) = \mathbf{x}_{\perp}^0. \end{array} \right. $$

Here, we are of course interested in the behavior of the sequence solution $(f^{\varepsilon})_{\varepsilon>0}$ to the Vlasov-Poisson system (1.3) when $\varepsilon \to 0$, which corresponds to the gyro-kinetic approximation of the
Vlasov-Poissons system. Following the work of L. Saint-Raymond [34], it can be proved — at least when $\mathbf{B}_{\text{ext}}$ is homogeneous — that the charge density $(\rho^{\varepsilon})_{\varepsilon>0}$ converges to the solution to the guiding center approximation

$$ (1.7) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho = 0, \\ -\Delta_{\mathbf{x}_{\perp}} \phi = \rho, \end{array} \right. $$

where the velocity $\mathbf{U}$ is

$$ \mathbf{U} = \frac{\mathbf{E} \wedge \mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2}, \quad \mathbf{E} = -\nabla_{\mathbf{x}_{\perp}} \phi. $$

We observe that the limit system (1.6) corresponds to the characteristic curves to the limit equation (1.7).

We seek a method that is able to capture these properties, while the numerical parameters may be kept independent of the stiffness degree of these scales. This concept is known and widely studied for dissipative systems in the framework of asymptotic preserving schemes [28, 29]. Contrary to collisional kinetic equations in hydrodynamic or diffusion asymptotic, collisionless equations like the Vlasov-Poisson system (1.3) involve time oscillations. In this context, the situation is more complicated than the one encountered in collisional regimes since we cannot expect any dissipative phenomenon. Therefore, the notion of two-scale convergence has been introduced both at the theoretical and numerical level [12, 21, 22] in order to derive asymptotic models. However, these asymptotic models, obtained after removing the fast scales, are valid only when $\varepsilon$ is small. We refer to E. Frénod and E. Sonnendrücker [22], and F. Golse and L. Saint-Raymond [24, 35] for a theoretical point of view on these questions, and E. Frénod, F. Salvarani and E. Sonnendrücker [21] for numerical applications of such techniques.

Another approach is to combine both disparate scales into one and single model. Such a decomposition can be done using a micro-macro approach (see [12] and the references therein). Such a model may be used when the small parameter of the equation is not everywhere small. Hence, a scheme for a micro-macro model can switch from one regime to another without any treatment of the transition between the different regimes. A different method consists in separating fast and slow time scales when such a structure can be identified [14] or [20].

These techniques work well when the magnetic field is uniform since fast scales can be computed from a formal asymptotic analysis, but for more complicated problems, that is, when the external magnetic field depends on time and position $\mathbf{x}$, the generalization of this approach is an open problem.

In this paper, we propose an alternative to such methods allowing to make direct simulations of systems (1.3) with large time steps with respect to $\varepsilon$. We develop numerical schemes that are able to deal with a wide range of values for $\varepsilon$, so-called Asymptotic Preserving (AP) class [29, 28], such schemes are consistent with the kinetic model for all positive value of $\varepsilon$, and degenerate into consistent schemes with the asymptotic model when $\varepsilon \to 0$.

Before presenting our time discretization technique, let us first briefly review the basic tools of particle-in-cell methods which are widely used for plasma physics simulations [5].

## 2. A BRIEF REVIEW OF PARTICLE METHODS

The numerical resolution of the Vlasov equation and related models is usually performed by Particle-In-Cell (PIC) methods which approximate the plasma by a finite number of particles. Trajectories of these particles are computed from characteristic curves (1.4) corresponding to the the Vlasov equation (1.3), whereas self-consistent fields are computed on a mesh of the physical space.
This method yields satisfying results with a relatively small number of particles but it is sometimes subject to fluctuations, due to the numerical noise, which are difficult to control. To improve the accuracy, direct numerical simulation techniques have been developed. The Vlasov equation is discretized in phase space using either semi-Lagrangian [17, 18, 32, 36], finite difference [39] or discontinuous Galerkin [4, 26] schemes. But these direct methods are very costly, hence several variants of particle methods have been developed over the past decades. In the Complex Particle Kinetic scheme introduced by Bateson and Hewett [2, 27], particles have a Gaussian shape that is transformed by the local shearing of the flow. Moreover they can be fragmented to probe for emerging features, and merged where fine particles are no longer needed. In the Cloud in Mesh (CM) scheme of Alard and Colombi [1] particles also have Gaussian shapes, and they are deformed by local linearization of the force field. More recently in [7], the authors proposed a Linearly-Transformed Particle-In-Cell method, that employs linear deformations of the particles.

Here we focus on the time discretization technique, hence we will only consider standard particle method even if our approach is completely independent on the choice of the particle method.

The particles method consists in approximating the initial condition $f_0$ in (1.3) by the following Dirac mass sum

$$f_N^0(\mathbf{x}, \mathbf{v}) := \sum_{1 \le k \le N} \omega_k \delta(\mathbf{x} - \mathbf{x}_k^0) \delta(\mathbf{v} - \mathbf{v}_k^0),$$

where $(\mathbf{x}_k^0, \mathbf{v}_k^0)_{1 \le k \le N}$ is a beam of $N$ particles distributed in the four dimensional phase space according to the density function $f_0$. Afterwards, one approximates the solution of (1.3), by

$$f_N(t, \mathbf{x}, \mathbf{v}) := \sum_{1 \le k \le N} \omega_k \delta(\mathbf{x} - \mathbf{X}_k(t)) \delta(\mathbf{v} - \mathbf{V}_k(t)),$$

where $(\mathbf{X}_k, \mathbf{V}_k)_{1 \le k \le N}$ is the position in phase space of particle $k$ moving along the characteristic curves (1.4) with the initial data $(\mathbf{x}_k^0, \mathbf{v}_k^0)$, for $1 \le k \le N$.

However when the Vlasov equation is coupled with the Poisson equation for the computation of the electric field, the Dirac mass has to be replaced by a smooth function $\varphi_\alpha$

$$f_{N,\alpha}^0(\mathbf{x}, \mathbf{v}) := \sum_{1 \le k \le N} \omega_k \varphi_\alpha(\mathbf{x} - \mathbf{x}_k^0) \varphi_\alpha(\mathbf{v} - \mathbf{v}_k^0),$$

where $\varphi_\varepsilon = \varepsilon^{-d}\varphi(\cdot/\varepsilon)$ is a particle shape function with radius proportional to $\varepsilon$, usually seen as a smooth approximation of the Dirac measure obtained by scaling a compactly supported “cut-off” function $\varphi$ for which common choices include B-splines and smoothing kernels with vanishing moments, see e.g. [30, 11].

Particle centers are then pushed forward at each time $t$ by following a numerical approximation of the flow (1.4), leading to

$$f_{N,\alpha}(t, \mathbf{x}, \mathbf{v}) := \sum_{1 \le k \le N} \omega_k \varphi_\alpha(\mathbf{x} - \mathbf{X}_k(t)) \varphi_\alpha(\mathbf{v} - \mathbf{V}_k(t)).$$

In the classical error analysis [3, 33], the above process is seen as

* An approximation (in the distribution sense) of the initial data by a collection of weighted Dirac measures;

* The exact transport of the Dirac particles along the flow;

* The smoothing of the resulting distribution

$$\sum_{1 \le k \le N} \omega_k \delta(\cdot - \mathbf{X}_k(t)) \delta(\cdot - \mathbf{V}_k(t))$$
with the convolution kernel $\varphi_{\varepsilon}$.

The classical error estimate reads then as follows [9]:

**Proposition 2.1.** Consider the Vlasov equation with a given electromagnetic field $(\mathbf{E}, \mathbf{B}_{\text{ext}})$ and a smooth initial datum $f^0 \in C_c^s(\mathbb{R}^d)$, with $s \ge 1$.

If for some prescribed integers $m > 0$ and $r > 0$, the cut-off $\varphi \ge 0$ has $m$-th order smoothness and satisfies a moment condition of order $r$, namely,

$$ \int_{\mathbb{R}^d} \varphi(\mathbf{y}) \, d\mathbf{y} = 1, \quad \int_{\mathbb{R}^d} |\mathbf{y}|^r \varphi(y) \, d\mathbf{y} < \infty, $$

and

$$ \int_{\mathbb{R}^d} y_1^{s_1} \dots y_d^{s_d} \varphi(\mathbf{y}) \, d\mathbf{y} = 0, \quad \text{for } \mathbf{s} = (s_1, \dots, s_d) \in \mathbb{N}^d \text{ with } 1 \le s_1 + \dots + s_d \le r - 1. $$

Then there exists a constant $C$ independent of $f_0$, $N$ or $\alpha$, such that we have for all $1 \le p \le +\infty$,

$$ \|f(t) - f_{N,\alpha}(t)\|_{L^p} \to 0, $$

when $N \to \infty$ and $\alpha \to 0$ where the ratio $N^{1/d}\alpha \ll 1$.

Note that following [9], it is also possible to get explicit order of convergence for the linear Vlasov equation. Let us also mention related papers where the convergence of a numerical scheme for the Vlasov-Poisson system is investigated. G.-H. Cottet and P.-A. Raviart [10] present a precise mathematical analysis of the particle method for solving the one-dimensional Vlasov-Poisson system. We also mention the papers of S. Wollman and E. Ozizmir [37] and S. Wollman [38] on the topic. K. Ganguly and H.D. Victory give a convergence result for the Vlasov-Maxwell system [23].

The rest of the paper is organized as follows. In Section 3 we present several time discretization techniques based on high-order semi-implicit schemes [6] for the Vlasov-Poisson system with a strong external magnetic field, and we prove uniform consistency of the schemes in the limit $\varepsilon \to 0$ with preservation of the order of accuracy (from first to third order accuracy). In Section 4 we perform a rigorous analysis of the first order scheme for smooth electromagnetic fields.

Section 5 is then devoted to numerical simulations for one single particle motion and for the Vlasov-Poisson model for various asymptotics $\varepsilon \approx 1$ and $\varepsilon \ll 1$, which illustrate the advantage of high order schemes.

### 3. A PARTICLE METHOD FOR VLASOV-POISSON SYSTEM WITH A STRONG MAGNETIC FIELD

Let us now consider the system (1.3) and apply a particle method, where the key issue is to design a uniformly stable scheme with respect to the parameter $\varepsilon > 0$, which is related to the magnitude of the external magnetic field. Assume that at time $t^n = n \Delta t$, the set of particles are located in $(\mathbf{x}_k^n, \mathbf{v}_k^n)_{1 \le k \le N}$, we want to solve the following system on the time interval $[t^n, t^{n+1}]$,

$$ (3.1) \qquad \left\{ \begin{array}{l} \varepsilon \frac{d\mathbf{X}_k}{dt} = \mathbf{V}_k, \\ \\ \varepsilon \frac{d\mathbf{V}_k}{dt} = \frac{1}{\varepsilon} \mathbf{V}_k \wedge \mathbf{B}_{\text{ext}}(t, \mathbf{X}_k) + \mathbf{E}(t, \mathbf{X}_k), \\ \\ \mathbf{X}(t^n) = \mathbf{x}_k^n, \quad \mathbf{V}(t^n) = \mathbf{v}_k^n, \end{array} \right. $$

where the electric field is computed from a discretization of the Poisson equation (1.3) on a mesh of the physical space.
The numerical scheme that we describe here is proposed in the framework of Particle-In-Cell method, where the solution $f$ is discretized as follows

$$f_{N,\alpha}^{n+1}(\mathbf{x}, \mathbf{v}) := \sum_{1 \le k \le N} \omega_k \varphi_{\alpha}(\mathbf{x} - \mathbf{x}_k^{n+1}) \varphi_{\alpha}(\mathbf{v} - \mathbf{v}_k^{n+1}),$$

where $(\mathbf{x}_k^{n+1}, \mathbf{v}_k^{n+1})$ represents an approximation of the solution $\mathbf{X}_k(t^{n+1})$ and $\mathbf{V}_k(t^{n+1})$ to (3.1).

When the Vlasov equation (1.3) is coupled with the Poisson equation (1.2), the electric field is computed in a macro-particle position $\mathbf{x}_k^{n+1}$ at time $t^{n+1}$ as follows

* Compute the density $\rho$

$$\rho_{h,\varepsilon}^{n}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{Z}^d} w_{\mathbf{k}} \varphi_{\varepsilon}(\mathbf{x} - \mathbf{x}_{\mathbf{k}}^{n}).$$

* Solve a discrete approximation to (1.2)

$$-\Delta_h \phi^n(\mathbf{x}) = \rho_{h,\varepsilon}^n(\mathbf{x}).$$

* Interpolate the electric field with the same order of accuracy on the points $(\mathbf{x}_k^n)_{\mathbf{k} \in \mathbb{Z}^d}$.

To discretize the system (3.1), we apply the strategy developed in [6] based on semi-implicit solver for stiff problems. In the rest of this section, we propose several numerical schemes to the system (3.1) for which the index $k \in \{1, \dots, N\}$ will be omitted.

**3.1. A first order semi-implicit scheme.** We start with the simplest semi-implicit scheme for (3.1), which is a combination of backward and forward Euler scheme. It gives for a fixed time step $\Delta t > 0$ and a given electric field $\mathbf{E}$ and an external magnetic field $\mathbf{B}_{\text{ext}}$,

$$(3.2) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\mathbf{x}^{n+1} - \mathbf{x}^n}{\Delta t} = \frac{\mathbf{v}^{n+1}}{\varepsilon}, \\[1.5em] \displaystyle \frac{\mathbf{v}^{n+1} - \mathbf{v}^n}{\Delta t} = \frac{1}{\varepsilon} \left( \frac{\mathbf{v}^{n+1}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}^n) + \mathbf{E}(t^n, \mathbf{x}^n) \right). \end{array} \right.$$

Note that only the second equation on $\mathbf{v}^{n+1}$ is fully implicit and requires the inversion of a linear operator. Then, from $\mathbf{v}^{n+1}$ the first equation gives the value for the position $\mathbf{x}^{n+1}$.

**Proposition 3.1** (Consistency in the limit $\varepsilon \to 0$ for a fixed $\Delta t$). Let us consider a time step $\Delta t > 0$, a final time $T > 0$ and set $N_T = [T/\Delta t]$. Assume that the sequence $(\mathbf{x}_{\varepsilon}^n, \mathbf{v}_{\varepsilon}^n)_{0 \le n \le N_T}$ given by (3.2) is such that for all $1 \le n \le N_T$, $(\mathbf{x}_{\varepsilon}^n, \varepsilon\mathbf{v}_{\varepsilon}^n)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon > 0$ and $(\mathbf{x}_{\varepsilon}^0, \varepsilon\mathbf{v}_{\varepsilon}^0)_{\varepsilon>0}$ converges in the limit $\varepsilon \to 0$ to some $(\mathbf{y}^0, 0)$. Then, for $1 \le n \le N_T$, $\mathbf{x}_{\varepsilon}^n \to \mathbf{y}^n$, as $\varepsilon \to 0$ and the limit $(\mathbf{y}^n)_{1 \le n \le N_T}$ is a consistent first order approximation with respect to $\Delta t$ of the guiding center equation provided by the scheme

$$(3.3) \qquad \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^n, \mathbf{y}^n)}{\|\mathbf{B}\|^2}.$$

*Proof.* For all $1 \le n \le N_T$, we consider $(\mathbf{x}_{\varepsilon}^n, \mathbf{v}_{\varepsilon}^n)$ the solution to (3.2) now labeled with respect to $\varepsilon > 0$. Since, the sequence $(\mathbf{x}_{\varepsilon}^n)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon > 0$, we can extract a subsequence still abusively labeled by $\varepsilon$ and find some $(\mathbf{y}^n)$ such that $\mathbf{x}_{\varepsilon}^n \to \mathbf{y}^n$ as $\varepsilon$ goes to zero. Then, we observe that the second equation of (3.2) can be written as

$$\varepsilon^2 \frac{\mathbf{v}_{\varepsilon}^{n+1} - \mathbf{v}_{\varepsilon}^n}{\varepsilon \Delta t} = \left( \frac{\mathbf{v}_{\varepsilon}^{n+1}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}_{\varepsilon}^n) + \mathbf{E}(t^n, \mathbf{x}_{\varepsilon}^n) \right).$$
and that, for any $0 \le n \le N_T$, $(\varepsilon \mathbf{v}_{\varepsilon})_{\varepsilon>0}$ is uniformly bounded. From this we conclude first that for any $1 \le n \le N_T$, $(\varepsilon^{-1} \mathbf{v}_{\varepsilon}^n)_{\varepsilon>0}$ is uniformly bounded then that $\varepsilon \mathbf{v}_{\varepsilon}^n \to 0$ for any $0 \le n \le N_T$. Therefore, taking the limit $\varepsilon \to 0$, it yields that for $0 \le n \le N_T - 1$,

$$ \varepsilon^{-1} \mathbf{v}_{\varepsilon}^{n+1} \to \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^n, \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2}, \text{ when } \varepsilon \to 0. $$

Substituting the limit of $\varepsilon^{-1} \mathbf{v}^{n+1}$ in the first equation of (3.2) we prove that the limit $\mathbf{y}^n$ satisfies (3.3). Since the limit point $\mathbf{y}_n$ is uniquely determined, actually all the sequence $(\mathbf{x}_{\varepsilon}^n)_{\varepsilon>0}$ converges. $\square$

**Remark 3.2.** The consistency provided by the latter result is far from being uniform with respect to the time step. However we do prove in the next section that the solution to (3.2) is both uniformly stable and consistent with respect to $\Delta t$ and $\varepsilon > 0$.

Of course, such a first order scheme is not accurate enough to describe correctly the long time behavior of the numerical solution, but it has the advantage of the simplicity and we will prove in the next section that it is uniformly stable with respect to the parameter $\varepsilon$ and the sequence $(\mathbf{x}_{\varepsilon}^n)$ converges to a consistent approximation of the guiding center model when $\varepsilon \to 0$.

Now, let us see how to generalize such an approach to second and third order schemes.

**3.2. Second order semi-implicit Runge-Kutta schemes.** Now, we consider second order schemes with two stages.

**3.2.1. A second order A-stable scheme.** A first example of scheme satisfying the second order conditions is given by a combination of Heun method (explicit part) and an A-stable second order singly diagonal implicit Runge-Kutta SDIRK method (implicit part) [25, 6]. The first stage corresponds to

$$ (3.4) \qquad \left\{ \begin{aligned} \mathbf{x}^{(1)} &= \mathbf{x}^n + \frac{\Delta t}{2\varepsilon} \mathbf{v}^{(1)}, \\ \mathbf{v}^{(1)} &= \mathbf{v}^n + \frac{\Delta t}{2\varepsilon} \left[ \frac{\mathbf{v}^{(1)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}^n) + \mathbf{E}(t^n, \mathbf{x}^n) \right]. \end{aligned} \right. $$

Then the second stage is given by

$$ (3.5) \qquad \left\{ \begin{aligned} \mathbf{x}^{(2)} &= \mathbf{x}^n + \frac{\Delta t}{2\varepsilon} \mathbf{v}^{(2)}, \\ \mathbf{v}^{(2)} &= \mathbf{v}^n + \frac{\Delta t}{2\varepsilon} \left[ \frac{\mathbf{v}^{(2)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^{n+1}, 2\mathbf{x}^{(1)} - \mathbf{x}^n) + \mathbf{E}(t^{n+1}, 2\mathbf{x}^{(1)} - \mathbf{x}^n) \right]. \end{aligned} \right. $$

Finally, the numerical solution at the new time step is

$$ (3.6) \qquad \left\{ \begin{aligned} \mathbf{x}^{n+1} &= \mathbf{x}^{(1)} + \mathbf{x}^{(2)} - \mathbf{x}^n, \\ \mathbf{v}^{n+1} &= \mathbf{v}^{(1)} + \mathbf{v}^{(2)} - \mathbf{v}^n. \end{aligned} \right. $$

A similar numerical scheme has been proposed in the framework of $\delta f$ simulation of the Vlasov-Poisson system [8].

Under stability assumptions on the numerical solution to (3.4)-(3.6), we get the following consistency result in the limit $\varepsilon \to 0$.
**Proposition 3.3** (Consistency in the limit $\epsilon \to 0$ for a fixed $\Delta t$). Let us consider a time step $\Delta t > 0$, a final time $T > 0$ and set $N_T = [T/\Delta t]$. Assume that the sequence $(\mathbf{x}_\epsilon^n, \mathbf{v}_\epsilon^n)_{0\le n\le N_T}$ given by (3.4)-(3.6) is such that for all $1 \le n \le N_T$, $(\mathbf{x}_\epsilon^n, \epsilon\mathbf{v}_\epsilon^n)_{\epsilon>0}$ is uniformly bounded with respect to $\epsilon > 0$ and $(\mathbf{x}_\epsilon^0, \epsilon\mathbf{v}_\epsilon^0)_{\epsilon>0}$ converges in the limit $\epsilon \to 0$ to some $(\mathbf{y}^0, 0)$. Then, for $1 \le n \le N_T$, $\mathbf{x}_\epsilon^n \to \mathbf{y}^n$, as $\epsilon \to 0$ and the limit $(\mathbf{y}^n)_{n\ge 1}$ is a consistent and second order approximation with respect to $\Delta t$ of the guiding center equation given by the scheme

$$ (3.7) \qquad \left\{ \begin{aligned} \mathbf{y}^{(1)} &= \mathbf{y}^n + \frac{\Delta t}{2} \left( \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^n, \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2} \right), \\ \mathbf{y}^{(2)} &= \mathbf{y}^n + \frac{\Delta t}{2} \left( \mathbf{E}(t^{n+1}, 2\mathbf{y}^{(1)} - \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^{n+1}, 2\mathbf{y}^{(1)} - \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2} \right), \end{aligned} \right. $$

and $\mathbf{y}^{n+1} = \mathbf{y}^{(1)} + \mathbf{y}^{(2)} - \mathbf{y}^n$.

*Proof.* We follow the lines of the proof of Proposition 3.1 and first choose a subsequence of $(\mathbf{x}_\epsilon^n)$ converging to some $\mathbf{y}^n$. The second equation of (3.4) implies that $\epsilon^{-1}\mathbf{v}_\epsilon^{(1)}$ is bounded. From the first equation of (3.4) it follows that so is $\mathbf{x}_\epsilon^{(1)}$. Then the second equation of (3.5) yields boundedness of $\epsilon^{-1}\mathbf{v}_\epsilon^{(2)}$ and the first that $\mathbf{x}_\epsilon^{(2)}$ is also bounded. Now from the second equation of (3.6) stems that $\epsilon\mathbf{v}_\epsilon^n$ converges to zero as $\epsilon \to 0$ for any $0 \le n \le N_T$. Coming back to the second equation of (3.4) we show that for all $n$, $\epsilon^{-1}\mathbf{v}_\epsilon^{(1)}$ converges when $\epsilon \to 0$ and

$$ \frac{\mathbf{v}_\epsilon^{(1)}}{\epsilon} \to \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^n, \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2}, \quad \text{when } \epsilon \to 0. $$

Using the first equation of (3.4) we conclude that $\mathbf{x}^{(1)}$ also converges and that its limit $\mathbf{y}^{(1)}$ is given by the first equation of (3.7). Going on with the same arguments shows that

$$ \frac{\mathbf{v}_\epsilon^{(2)}}{\epsilon} \to \mathbf{E}(t^{n+1}, 2\mathbf{y}^{(1)} - \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^{n+1}, 2\mathbf{y}^{(1)} - \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2}, \quad \text{when } \epsilon \to 0 $$

and that $\mathbf{x}^{(2)}$ converges to a $\mathbf{y}^{(2)}$ given by the second equation of (3.7). One may then take a limit in the first equation of (3.6) and conclude that indeed the limit $\mathbf{y}^n$ satisfies (3.7). Again uniqueness supplies the convergence of the whole sequence. $\square$

### 3.2.2. A second order L-stable scheme.
Another choice is a combination of Runge-Kutta method (explicit part) and an L-stable second order SDIRK method in the implicit part. This implicit scheme should give better stability properties on the numerical solution with respect to the stiffness parameter $\epsilon > 0$.

We first choose $\gamma > 0$ as the smallest root of the polynomial $\gamma^2 - 2\gamma + 1/2 = 0$, i.e. $\gamma = 1 - 1/\sqrt{2}$, then the scheme is given by the following two stages. First, we have

$$ (3.8) \qquad \left\{ \begin{aligned} \mathbf{x}^{(1)} &= \mathbf{x}^n + \frac{\gamma \Delta t}{\epsilon} \mathbf{v}^{(1)}, \\ \mathbf{v}^{(1)} &= \mathbf{v}^n + \frac{\gamma \Delta t}{\epsilon} \mathbf{F}^{(1)}, \end{aligned} \right. $$

with

$$ \mathbf{F}^{(1)} := \frac{\mathbf{v}^{(1)}}{\epsilon} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}^n) + \mathbf{E}(t^n, \mathbf{x}^n). $$

For the second stage, we first define

$$ (3.9) \qquad \hat{t}^{(1)} := t^n + \frac{\Delta t}{2\gamma}, \quad \hat{\mathbf{x}}^{(1)} := \mathbf{x}^n + \frac{\Delta t}{2\gamma\epsilon}\mathbf{v}^{(1)}, $$
then the solution $(\mathbf{x}^{(2)}, \mathbf{v}^{(2)})$ is given by

$$ (3.10) \qquad \left\{ \begin{array}{l} \mathbf{x}^{(2)} = \mathbf{x}^n + \frac{(1-\gamma)\Delta t}{\varepsilon} \mathbf{v}^{(1)} + \frac{\gamma\Delta t}{\varepsilon} \mathbf{v}^{(2)}, \\ \mathbf{v}^{(2)} = \mathbf{v}^n + \frac{(1-\gamma)\Delta t}{\varepsilon} \mathbf{F}^{(1)} + \frac{\gamma\Delta t}{\varepsilon} \mathbf{F}^{(2)}, \end{array} \right. $$

with

$$ \mathbf{F}^{(2)} := \frac{\mathbf{v}^{(2)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(\hat{t}^{(1)}, \hat{\mathbf{x}}^{(1)}) + \mathbf{E}(\hat{t}^{(1)}, \hat{\mathbf{x}}^{(1)}). $$

Finally, the numerical solution at the new time step is

$$ (3.11) \qquad \left\{ \begin{array}{l} \mathbf{x}^{n+1} = \mathbf{x}^{(2)}, \\ \mathbf{v}^{n+1} = \mathbf{v}^{(2)}. \end{array} \right. $$

Under stability assumptions on the numerical solution to (3.4)-(3.6), we get the following consistency result in the limit $\varepsilon \to 0$.

**Proposition 3.4** (Consistency in the limit $\varepsilon \to 0$ for a fixed $\Delta t$). Let us consider a time step $\Delta t > 0$, a final time $T > 0$ and set $N_T = [T/\Delta t]$. Assume that the sequence $(\mathbf{x}_\varepsilon^n, \mathbf{v}_\varepsilon^n)_{0 \le n \le N_T}$ given by (3.8)-(3.11) is such that for all $1 \le n \le N_T$, $(\mathbf{x}_\varepsilon^n, \varepsilon\mathbf{v}_\varepsilon^n)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon > 0$ and $(\mathbf{x}_\varepsilon^0, \varepsilon\mathbf{v}_\varepsilon^0)_{\varepsilon>0}$ converges in the limit $\varepsilon \to 0$ to some $(\mathbf{y}^0, 0)$. Then, for $1 \le n \le N_T$, $\mathbf{x}_\varepsilon^n \to \mathbf{y}^n$, as $\varepsilon \to 0$ and the limit $(\mathbf{y}^n)_{n \ge 1}$ is a consistent second order approximation of the guiding center equation, given by

$$ (3.12) \qquad \left\{ \begin{array}{l} \mathbf{U}^n = \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^n, \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2}, \\ \mathbf{y}^{n+1} = \mathbf{y}^n + (1-\gamma)\Delta t \mathbf{U}^n + \gamma\Delta t \mathbf{U}^{(1)}, \end{array} \right. $$

where

$$ \hat{\mathbf{y}}^{(1)} := \mathbf{y}^n + \frac{\Delta t}{2\gamma} \mathbf{U}^n, \quad \mathbf{U}^{(1)} := \mathbf{E}(\hat{t}^{(1)}, \hat{\mathbf{y}}^{(1)}) \wedge \frac{\mathbf{B}_{\text{ext}}(\hat{t}^{(1)}, \hat{\mathbf{y}}^{(1)})}{\|\mathbf{B}_{\text{ext}}\|^2}. $$

We omit the proof of Proposition 3.4 as almost identical to the one of Proposition 3.3.

The present scheme is L-stable, which means uniformly linearly stable with respect to $\Delta t$.

**3.3. Third order semi-implicit Runge-Kutta schemes.** A third order semi-implicit scheme is given by a four stages Runge-Kutta method introduced in the framework of hyperbolic systems with stiff source terms [6]. First, we set $\alpha = 0.24169426078821$, $\beta = \alpha/4$ and $\eta = 0.12915286960590$ and $\gamma = 1/2 - \alpha - \beta - \eta$. Then we construct the first stage as

$$ (3.13) \qquad \left\{ \begin{array}{l} \mathbf{x}^{(1)} = \mathbf{x}^n + \frac{\alpha \Delta t}{\varepsilon} \mathbf{v}^{(1)}, \\ \mathbf{v}^{(1)} = \mathbf{v}^n + \frac{\alpha \Delta t}{\varepsilon} \mathbf{F}^{(1)}, \end{array} \right. $$

with

$$ \mathbf{F}^{(1)} := \frac{\mathbf{v}^{(1)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}^n) + \mathbf{E}(t^n, \mathbf{x}^n). $$
For the second stage, we have

$$ (3.14) \qquad \left\{ \begin{array}{l} \mathbf{x}^{(2)} = \mathbf{x}^n - \frac{\alpha \Delta t}{\varepsilon} \mathbf{v}^{(1)} + \frac{\alpha \Delta t}{\varepsilon} \mathbf{v}^{(2)}, \\ \mathbf{v}^{(2)} = \mathbf{v}^n - \frac{\alpha \Delta t}{\varepsilon} \mathbf{F}^{(1)} + \frac{\alpha \Delta t}{\varepsilon} \mathbf{F}^{(2)}, \end{array} \right. $$

with

$$ \mathbf{F}^{(2)} := \frac{\mathbf{v}^{(2)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}^n) + \mathbf{E}(t^n, \mathbf{x}^n). $$

Then, for the third stage we set

$$ (3.15) \qquad \left\{ \begin{array}{l} \mathbf{x}^{(3)} = \mathbf{x}^n + \frac{(1-\alpha)\Delta t}{\varepsilon} \mathbf{v}^{(2)} + \frac{\alpha \Delta t}{\varepsilon} \mathbf{v}^{(3)}, \\ \mathbf{v}^{(3)} = \mathbf{v}^n + \frac{(1-\alpha)\Delta t}{\varepsilon} \mathbf{F}^{(2)} + \frac{\alpha \Delta t}{\varepsilon} \mathbf{F}^{(3)}, \end{array} \right. $$

with

$$ \left\{ \begin{array}{l} \mathbf{F}^{(3)} := \frac{\mathbf{v}^{(3)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^{n+1}, \bar{\mathbf{x}}^{(2)}) + \mathbf{E}(t^{n+1}, \bar{\mathbf{x}}^{(2)}), \\ \bar{\mathbf{x}}^{(2)} := \mathbf{x}^n + \frac{\Delta t}{\varepsilon} \mathbf{v}^{(2)}. \end{array} \right. $$

Finally, for the fourth stage we set

$$ (3.16) \qquad \left\{ \begin{array}{l} \mathbf{x}^{(4)} = \mathbf{x}^n + \frac{\beta \Delta t}{\varepsilon} \mathbf{v}^{(1)} + \frac{\eta \Delta t}{\varepsilon} \mathbf{v}^{(2)} + \frac{\gamma \Delta t}{\varepsilon} \mathbf{v}^{(3)} + \frac{\alpha \Delta t}{\varepsilon} \mathbf{v}^{(4)}, \\ \mathbf{v}^{(4)} = \mathbf{v}^n + \frac{\beta \Delta t}{\varepsilon} \mathbf{F}^{(1)} + \frac{\eta \Delta t}{\varepsilon} \mathbf{F}^{(2)} + \frac{\gamma \Delta t}{\varepsilon} \mathbf{F}^{(3)} + \frac{\alpha \Delta t}{\varepsilon} \mathbf{F}^{(4)}, \end{array} \right. $$

with

$$ \left\{ \begin{array}{l} \mathbf{F}^{(4)} := \frac{\mathbf{v}^{(4)}}{\varepsilon} \wedge \mathbf{B}_{\text{ext}}(t^{n+1/2}, \bar{\mathbf{x}}^{(3)}) + \mathbf{E}(t^{n+1/2}, \bar{\mathbf{x}}^{(3)}), \\ \bar{\mathbf{x}}^{(3)} := \mathbf{x}^n + \frac{\Delta t}{4\varepsilon} (\mathbf{v}^{(2)} + \mathbf{v}^{(3)}), \end{array} \right. $$

and the numerical solution at the new time step is

$$ (3.17) \qquad \left\{ \begin{array}{l} \mathbf{x}^{n+1} = \mathbf{x}^n + \frac{\Delta t}{6\varepsilon} (\mathbf{v}^{(2)} + \mathbf{v}^{(3)} + 4\mathbf{v}^{(4)}), \\ \mathbf{v}^{n+1} = \mathbf{v}^n + \frac{\Delta t}{6\varepsilon} (\mathbf{F}^{(2)} + \mathbf{F}^{(3)} + 4\mathbf{F}^{(4)}). \end{array} \right. $$

As for the previous schemes, under uniform stability assumptions with respect to $\varepsilon > 0$, we prove the following Proposition

**Proposition 3.5** (Consistency in the limit $\varepsilon \to 0$ for a fixed $\Delta t$). Let us consider a time step $\Delta t > 0$, a final time $T > 0$ and set $N_T = [T/\Delta t]$. Assume that the sequence $(\mathbf{x}_\varepsilon^n, \mathbf{v}_\varepsilon^n)_{0 \le n \le N_T}$ given by (3.13)-(3.17) is such that for all $1 \le n \le N_T$, $(\mathbf{x}_\varepsilon^n, \varepsilon\mathbf{v}_\varepsilon^n)_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon > 0$ and $(\mathbf{x}_\varepsilon^0, \varepsilon\mathbf{v}_\varepsilon^0)_{\varepsilon>0}$ converges in the limit $\varepsilon \to 0$ to some $(\mathbf{y}^0, 0)$. Then, for $1 \le n \le N_T$, $\mathbf{x}_\varepsilon^n \to \mathbf{y}^n$, as $\varepsilon \to 0$ and the limit $(\mathbf{y}^n)_{n \ge 1}$ is a consistent third order approximation
of the guiding center equation provided by the scheme

$$ (3.18) \qquad \left\{ \begin{array}{l} \mathbf{U}^{(2)} = \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}(t^n, \mathbf{y}^n)}{\|\mathbf{B}_{\text{ext}}\|^2}, \\[1em] \mathbf{U}^{(3)} = \mathbf{E}(t^{n+1}, \mathbf{y}^{(2)}) \wedge \frac{\mathbf{B}_{\text{ext}}(t^{n+1}, \mathbf{y}^{(2)})}{\|\mathbf{B}_{\text{ext}}\|^2}, \\[1em] \mathbf{U}^{(4)} = \mathbf{E}(t^{n+1/2}, \mathbf{y}^{(3)}) \wedge \frac{\mathbf{B}_{\text{ext}}(t^{n+1/2}, \mathbf{y}^{(3)})}{\|\mathbf{B}_{\text{ext}}\|^2}, \end{array} \right. $$

where

$$ \left\{ \begin{array}{l} \mathbf{y}^{(2)} = \mathbf{y}^n + \Delta t \mathbf{U}^{(2)}, \\ \mathbf{y}^{(3)} = \mathbf{y}^n + \frac{\Delta t}{4} (\mathbf{U}^{(2)} + \mathbf{U}^{(3)}), \end{array} \right. $$

and

$$ \mathbf{y}^{n+1} = \mathbf{y}^n + \frac{\Delta t}{6} (\mathbf{U}^{(2)} + \mathbf{U}^{(3)} + 4 \mathbf{U}^{(4)}). $$

*Proof.* We follow the lines of the proof of Proposition 3.3. In particular we prove successively the boundedness of $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(1)}$, $\mathbf{x}_{\varepsilon}^{(1)}$, $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(2)}$, $\mathbf{x}_{\varepsilon}^{(2)}$, $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(3)}$, $\mathbf{x}_{\varepsilon}^{(3)}$, $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(4)}$, $\mathbf{x}_{\varepsilon}^{(4)}$. Arguing inductively on $n$ one then shows that $\varepsilon\mathbf{v}_{\varepsilon}^{n}$, $\mathbf{F}^{(1)}$, $\mathbf{F}^{(2)}$, $\mathbf{F}^{(3)}$, $\mathbf{F}^{(4)}$ all converge to zero when $\varepsilon \to 0$. Then for any converging subsequence of $\mathbf{x}_{\varepsilon}^{n}$ we may successively identify limits for $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(1)}$, $\mathbf{x}_{\varepsilon}^{(1)}$, $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(2)}$, $\mathbf{x}_{\varepsilon}^{(2)}$, $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(3)}$, $\mathbf{x}_{\varepsilon}^{(3)}$, $\varepsilon^{-1}\mathbf{v}_{\varepsilon}^{(4)}$, $\mathbf{x}_{\varepsilon}^{(4)}$ and prove that any accumulating point of $\mathbf{x}_{\varepsilon}^{n}$ solves the limiting scheme involving (3.18). Hence the result. $\square$

## 4. ANALYSIS OF THE FIRST-ORDER SEMI-IMPLICIT SCHEME

Consider the first order Euler semi-implicit scheme

$$ (4.1) \qquad \left\{ \begin{array}{l} \varepsilon \frac{\mathbf{x}_k^{n+1} - \mathbf{x}_k^n}{\Delta t} = \mathbf{v}_k^{n+1}, \\[1em] \varepsilon \frac{\mathbf{v}_k^{n+1} - \mathbf{v}_k^n}{\Delta t} = \frac{1}{\varepsilon} \mathbf{v}_k^{n+1} \wedge \mathbf{B}_{\text{ext}}(t^n, \mathbf{x}_k^n) + \mathbf{E}(t^n, \mathbf{x}_k^n), \\[1em] \mathbf{x}_k^0 = \mathbf{x}_k^0, \quad \mathbf{v}_k^0 = \mathbf{v}_k^0. \end{array} \right. $$

We focus here on the case where both $\mathbf{E}$ and $\mathbf{B}_{\text{ext}}$ are given external fields. Since, in this case, there is no coupling between particles we may focus safely on one of them and drop the $k$ suffix. Then we prove

**Theorem 4.1** (Uniform consistency with respect to $\varepsilon$). We set

$$ \lambda := \frac{\Delta t}{\varepsilon^2}, $$

assume that the $(\mathbf{E}, \mathbf{B}_{\text{ext}}) \in W_{\text{loc}}^{1,\infty}(\mathbb{R}^+ \times \mathbb{R}^2)$ and consider the solution $(\mathbf{x}^n, \mathbf{v}^n)$ to (4.1). Then there exist positive constants $C$ and $\lambda_0$, only depending on $\mathbf{E}$, $\mathbf{B}_{\text{ext}}$ and $t^n$ such that when $1 \le \Delta t \le 1$ and $\lambda \ge \lambda_0$

$$ \|\mathbf{x}^n - \mathbf{y}^n\| \le \frac{C \Delta t}{\lambda} [1 + \|\varepsilon^{-1}\mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)]\|], $$

<footnote n="1">¹Obviously one may replace 1 by any upper bound on $\Delta t$ but $\lambda_0$ and $C$ would depend on this upper bound.</footnote>
where $(\mathbf{y}^n)$ corresponds to the numerical solution of the guiding center model

$$ (4.2) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = \mathbf{E}(t^n, \mathbf{y}^n) \wedge \frac{\mathbf{B}_{\text{ext}}}{\|\mathbf{B}_{\text{ext}}\|^2}, \\ \\ \mathbf{y}^0 = \mathbf{x}^0. \end{array} \right. $$

*Proof.* To start with we analyze the case where $\mathbf{E}$ is a given bounded Lipschitz field and $\mathbf{B}_{\text{ext}}$ is constant and of norm 1. Our estimates shall be expressed in terms of

$$ K_0 = \|\mathbf{E}\|_{L^\infty} \qquad K_t = \left\| \frac{\partial \mathbf{E}}{\partial t} \right\|_{L^\infty} \qquad K_x = \|\nabla_x \mathbf{E}\|_{L^\infty}. $$

For conclusion's sake we also introduce the operator to denote $R[\mathbf{W}] = \mathbf{W} \wedge \mathbf{B}_{\text{ext}}$.

By introducing the key quantity for $n \ge 1$,

$$ \mathbf{z}^n = \varepsilon^{-1} \mathbf{v}^n - R[\mathbf{E}(t^{n-1}, \mathbf{x}^{n-1})], $$

the first equation of (4.1) reads

$$ \mathbf{x}^n = \mathbf{x}^{n-1} + \Delta t R[\mathbf{E}(t^{n-1}, \mathbf{x}^{n-1})] + \Delta t \mathbf{z}^n, \quad n \ge 1 $$

while $\mathbf{z}^1 = [\text{Id} - \lambda R]^{-1}(\varepsilon^{-1} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)])$ and

$$ \mathbf{z}^{n+1} = [\text{Id} - \lambda R]^{-1}(\mathbf{z}^n - R[\mathbf{E}(t^n, \mathbf{x}^n) - \mathbf{E}(t^{n-1}, \mathbf{x}^{n-1})]), \quad n \ge 1. $$

We first observe that $\|R\| \le 1$ and $[\text{Id} - \lambda R]^{-1} = (\text{Id} + \lambda R)/(1 + \lambda^2)$ since $R^2 = -\text{Id}$.

This leads to

$$ \|\mathbf{z}^{n+1}\| \le \frac{1+\lambda}{1+\lambda^2} \left[ (1+K_x \Delta t) \|\mathbf{z}^n\| + (K_t + K_x K_0) \Delta t \right], \quad n \ge 1. $$

Assuming

$$ \frac{1+\lambda}{1+\lambda^2} (1+K_x \Delta t) < 1 $$

and introducing

$$ a = \frac{1+\lambda}{1+\lambda^2} \frac{K_t + K_x K_0}{1 - \frac{1+\lambda}{1+\lambda^2}(1+K_x \Delta t)} \Delta t $$

we infer

$$ \|\mathbf{z}^n\| \le a + \left( \frac{1+\lambda}{1+\lambda^2} (1+K_x \Delta t) \right)^{n-1} b, \quad n \ge 1, $$

where

$$ b = a + \frac{1+\lambda}{1+\lambda^2} \left\| \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right\|. $$

For comparison we define $\mathbf{y}^n$ solution to (4.2). Then,

$$ \|\mathbf{x}^n - \mathbf{y}^n\| \le (1 + K_x \Delta t) \|\mathbf{x}^{n-1} - \mathbf{y}^{n-1}\| \\ \qquad + \Delta t \left[ a + \left( \frac{1+\lambda}{1+\lambda^2} (1+K_x \Delta t) \right)^{n-1} b \right], \quad n \ge 1. $$

Hence assuming moreover $K_x > 0$ (or replacing $K_x$ with some arbitrary positive number if $\|\nabla_x \mathbf{E}\|_{L^\infty} = 0$)

$$ \|\mathbf{x}^n - \mathbf{y}^n\| \le \frac{a}{K_x} [1 + (1 + K_x \Delta t)^n] \\ \qquad + \frac{b \Delta t}{(1 + K_x \Delta t)[1 - \frac{1+\lambda}{1+\lambda^2}]} \left[1 + \left(\frac{1+\lambda}{1+\lambda^2}\right)^n\right] (1 + K_x \Delta t)^n. $$
As a conclusion, for any $0 \le \theta < 1$, there exists $C_{\theta}$ such that if

$$ \frac{1 + \lambda}{1 + \lambda^2} (1 + K_x \Delta t) \le \theta $$

then, it yields

$$ (4.3) \qquad \|\mathbf{x}^n - \mathbf{y}^n\| \le C_{\theta} \frac{\Delta t}{\lambda} \left[ K_0 + \frac{K_t}{K_x} + \left\| \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right\| \right] e^{K_x n \Delta t}, $$

which concludes the proof of Theorem 4.1 when the magnetic field is uniform.

We now relax the assumption that $\mathbf{B}_{\text{ext}}$ is constant and of norm 1. We set $\mathbf{b}_{\text{ext}} = \|\mathbf{B}_{\text{ext}}\|$ and let $R$ be dependent on $t$ and $\mathbf{x}$. We observe that now $R^2 = -\mathbf{b}_{\text{ext}}^2 \text{Id}$ so that $[\text{Id} - \lambda R]^{-1} = (\text{Id} + \lambda R)/(1 + \lambda^2 \mathbf{b}_{\text{ext}}^2)$. Then introducing the drift force

$$ \mathbf{F}(t, \mathbf{x}) = \frac{1}{\|\mathbf{B}_{\text{ext}}(t, \mathbf{x})\|^2} \mathbf{E}(t, \mathbf{x}) \wedge \mathbf{B}_{\text{ext}}(t, \mathbf{x}) $$

we essentially obtain the same estimates with $\frac{1+\lambda}{1+\lambda^2}$ replaced by

$$ \frac{1}{\Lambda} := \left\| \frac{1 + \lambda \mathbf{b}_{\text{ext}}}{1 + \lambda^2 \mathbf{b}_{\text{ext}}^2} \right\|_{L^{\infty}} $$

and

$$ K_0 = \|\mathbf{F}\|_{L^{\infty}} \qquad K_t = \left\| \frac{\partial \mathbf{F}}{\partial t} \right\|_{L^{\infty}} \qquad K_x = \|\nabla_{\mathbf{x}} \mathbf{F}\|_{L^{\infty}}. $$

Indeed introducing for $n \ge 1$,

$$ \mathbf{z}^n = \frac{\mathbf{v}^n}{\varepsilon} - \mathbf{F}(t^{n-1}, \mathbf{x}^{n-1}), $$

the scheme is written

$$ \mathbf{z}^1 = [\text{Id} - \lambda R(t^0, \mathbf{x}^0)]^{-1} \left( \frac{\mathbf{v}^0}{\varepsilon} - \mathbf{F}(t^0, \mathbf{x}^0) \right), $$

and then

$$ \mathbf{z}^{n+1} = [\text{Id} - \lambda R(t^n, \mathbf{x}^n)]^{-1} (\mathbf{z}^n - (\mathbf{F}(t^n, \mathbf{x}^n) - \mathbf{F}(t^{n-1}, \mathbf{x}^{n-1}))), \quad n \ge 1 $$

together with

$$ \mathbf{x}^n = \mathbf{x}^{n-1} + \Delta t \mathbf{F}(t^{n-1}, \mathbf{x}^{n-1}) + \Delta t \mathbf{z}^n, \quad n \ge 1. $$

Mark that if $\lambda \times (\inf \|\mathbf{B}_{\text{ext}}\|) > \sqrt{2} - 1$ then

$$ \Lambda \ge \frac{1 + \lambda^2 \times (\inf \|\mathbf{B}_{\text{ext}}\|)^2}{1 + \lambda \times (\inf \|\mathbf{B}_{\text{ext}}\|)} . $$

□

Note that this result can be slightly improved when we modify the initial condition of the asymptotic discrete model. Indeed, consider $\mathbf{y}^n$ solving

$$ (4.4) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = \mathbf{E}(t^n, \mathbf{y}^n) \wedge \mathbf{B}_{\text{ext}}, \\ \mathbf{y}^0 = \mathbf{x}^0 + \varepsilon (\mathbf{v}^0 \wedge \mathbf{B}_{\text{ext}} + \varepsilon \mathbf{E}(t^0, \mathbf{x}^0)). \end{array} \right. $$

The gain is that now

$$ \|\mathbf{x}^1 - \mathbf{y}^1\| \le \frac{\Delta t}{\lambda} \left[ K_x \Delta t + \frac{1 + \lambda}{1 + \lambda^2} \right] \left\| \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right\| $$
since

$$ \mathbf{z}^1 = \frac{1}{\lambda} R \left[ \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right] - \frac{1}{\lambda} R [\mathrm{Id} - \lambda R]^{-1} \left[ \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right]. $$

This leads, for $n \ge 1$ to

$$ \begin{aligned} \|\mathbf{x}^n - \mathbf{y}^n\| &\le \frac{a}{K_x} \left[ 1 + (1 + K_x \Delta t)^{n-1} \right] \\ &\quad + \frac{b \Delta t}{1 - \frac{1+\lambda}{1+\lambda^2}} \frac{1+\lambda}{1+\lambda^2} \left[ 1 + \left( \frac{1+\lambda}{1+\lambda^2} \right)^{n-1} \right] (1 + K_x \Delta t)^{n-1} \\ &\quad + \frac{\Delta t}{\lambda} \left[ K_x \Delta t + \frac{1+\lambda}{1+\lambda^2} \right] \left\| \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right\| (1 + K_x \Delta t)^{n-1} \end{aligned} $$

then with notation as above there exists a constant $C = C_\theta (K_0 + \frac{K_t}{K_x} + K_x + 1)e^{K_x t^n} > 0$ such that

$$ (4.5) \qquad \|\mathbf{x}^n - \mathbf{y}^n\| \le C \frac{\Delta t}{\lambda} \left[ 1 + \left(\frac{1}{\lambda} + \Delta t\right) \left\| \frac{1}{\varepsilon} \mathbf{v}^0 - R[\mathbf{E}(t^0, \mathbf{x}^0)] \right\| \right]. $$

Concerning the analysis of high-order schemes, it is not straightforward to adapt directly the strategy of Theorem 4.1. Indeed, the use of a semi-implicit scheme does not necessarily guarantee that the particle trajectories are under control.

**Remark 4.2.** Consider the scheme (3.4)-(3.6) in the simplest situation where the electric field is zero and the magnetic field is $\mathbf{B}_{\mathrm{ext}} = (0, 0, 1)$, we show that the scheme preserves the kinetic energy and we have

$$ \|\mathbf{v}^n\|^2 = \|\mathbf{v}^0\|^2, \quad \forall n \in \mathbb{N}, $$

hence as $\varepsilon$ goes to zero, the velocity $\mathbf{v}^n/\varepsilon$ cannot converge to the null guiding center velocity except if the initial velocity does. Fortunately, for the other high order schemes (3.8)-(3.11) and (3.13)-(3.17), kinetic energy is dissipated and converges to 0.

## 5. NUMERICAL SIMULATIONS

In this section, we discuss some examples to validate and to compare the different time discretization schemes. We first consider the single motion of a particle under the effect of a given electromagnetic field. It allows us to illustrate the ability of the semi-implicit schemes to capture the guiding center velocity with large time step $\Delta t$ in the limit $\varepsilon \to 0$.

Then we consider the Vlasov-Poisson system with an external magnetic field. A classical particle-in-cell method is applied with different time discretization techniques to compute the particle trajectories. Hence this collection of charged particles move and give rise to a self-consistent electric field, obtained by solving numerically the Poisson equation in a space grid.

### 5.1. One single particle motion.
Before going to the statistical descriptions, let us investigate the accuracy and stability properties of the semi-implicit algorithms presented in Section 3 on the motion of individual particles in a given electromagnetic field.

Here we consider an electric field $\mathbf{E} = -\nabla\phi$, where

$$ \phi(\mathbf{x}) = \frac{1}{2} (\|\mathbf{x}\|^2 + \alpha \cos^2(2\pi y)), \quad \mathbf{x} = (x, y) \in \mathbb{R}^2, $$

with $\alpha = 0.02$ and a magnetic field $\mathbf{B}(\mathbf{x}) = 1 + 10^{-1} \sin(2\pi x)$ with $\mathbf{x} = (x, y) \in \mathbb{R}^2$. We choose for all simulations $\Delta t = 0.1$ and the initial data as $\mathbf{x}^0 = (1, 1.4)$ and $\varepsilon^{-1} \mathbf{v}^0 = (3, 5)$, such that the initial data $\mathbf{z}^0 = \varepsilon^{-1} \mathbf{v}^0 + \mathbf{F}(0, \mathbf{x}^0)$ is bounded with respect to $\varepsilon > 0$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. One single particle motion. Numerical solution obtained with a large time step $\Delta t = 0.1$ with the 1st order scheme (3.2), the 2nd order schemes (3.4)-(3.6) and (3.8)-(3.11) and the third order scheme (3.13)-(3.17) for $\varepsilon = 1$ (top) and $\varepsilon = 5 \cdot 10^{-1}$ (bottom):</figcaption></figure>

(a) particle trajectory in physical space ($\mathbf{x}^n$)$_{0\le n\le N_T}$

(b) particle velocity ($\mathbf{v}^n$)$_{0\le n\le N_T}$.

Thus, we apply the schemes proposed in Section 3 to compute a numerical solution to (3.1) and present the particle trajectory and velocity in Figures 1 and 2. These results are compared with those obtained with a fourth order Runge-Kutta scheme using a small time step.

In Figure 1, we first investigate the case where $\varepsilon$ is of order one ($\varepsilon = 1$ and $0.5$), which is the non stiff regime. On the left column, we clearly observe that the space trajectory obtained from high order schemes agree very well with the reference trajectory, whereas after few time steps the first order scheme does not give satisfying results. On the right hand side, we present the evolution of the velocity at each time step and compare it with the reference velocity and the guiding center velocity. In this non stiff regime, the velocity obtained from high order schemes
<figure><img src="image_3.png" /><figcaption>FIGURE 2. One single particle motion. Numerical solution obtained with a large time step $\Delta t = 0.1$ with the 1st order scheme (3.2), the 2nd order schemes (3.4)-(3.6) and (3.8)-(3.11) and the third order scheme (3.13)-(3.17) for $\varepsilon = 10^{-1}$ (top) and $\varepsilon = 10^{-2}$ (bottom):</figcaption></figure>

(a) particle trajectory in physical space $(\mathbf{x}^n)_{0 \le n \le N_T}$

(b) particle velocity $(\mathbf{v}^n)_{0 \le n \le N_T}$

coincides with the reference velocity and the guiding center velocity is meaningless. Therefore, this first test illustrates the ability of high order schemes to describe accurately the particle motion in phase space when $\varepsilon \sim 1$.

In Figure 2, we now propose the numerical results when $\varepsilon \ll 1$, that is, $\varepsilon = 0.1$ and $\varepsilon = 0.01$, which corresponds to the high field regime. In that case, the space trajectory in the orthogonal plane to the magnetic field can be decomposed into a relatively slow motion due to the guiding center velocity

$$ \mathbf{F}(t, \mathbf{x}) = \frac{1}{\|\mathbf{B}_{\text{ext}}(t, \mathbf{x})\|^2} \mathbf{E}(t, \mathbf{x}) \wedge \mathbf{B}_{\text{ext}}(t, \mathbf{x}) $$
along the field line and a fast circular motion with a frequency of order $1/\epsilon$. Our aim here is to capture the slow motion using a fixed time step $\Delta t = 0.1$ independently of the value $\epsilon \ll 1$.

On the one hand, the space trajectory (left column of Figure 2) of the numerical solution remains stable for various $\epsilon > 0$ even if we do not solve the fast scales. Moreover, when $\epsilon \to 0$, the numerical solution approaches the correct trajectory and fast fluctuations are somehow filtered thanks to the implicit treatment of the velocity $\mathbf{v}^n$. As in the previous simulations, we observe a discrepancy between the first order scheme and other high order schemes for large time.

Furthermore, we focus on the velocity variable $(\mathbf{v}^n)_{0 \le n \le N_T}$ and compare its time evolution with the reference velocity and the guiding center velocity (right column of Figure 2). Since we use a large time step $\Delta t = 0.1$, we cannot expect to follow all the details of the velocity variable, but only the slow motion corresponding to the guiding center velocity. Now we clearly observe a different behavior of the numerical solutions. On the one hand, the scheme (3.4)-(3.6), which preserves the kinetic energy $\frac{1}{2}|\mathbf{v}^n|^2$ when $\mathbf{E} = 0$ and $\mathbf{B}_{\text{ext}} = (0, 0, 1)$ (see Remark 4.2), still oscillates with a large amplitude and the velocity does not coincide with the guiding center velocity. On the other hand, the schemes (3.2), (3.8)-(3.11) and (3.13)-(3.17) are more dissipative and after few time steps the velocity follows the line of the guiding center velocity (see the zoom on the right column of Figure 2). Hence the amplitude of oscillations in the physical space are diminishing and the particle follows the trajectory corresponding to the guiding center model (1.7).

As a conclusion, these elementary numerical simulations confirm the ability of the semi-implicit discretization to capture the slow motion corresponding to the guiding center model (1.7) uniformly with respect to $\epsilon \ll 1$ and the interest of high order time discretization for the long time behavior of the solution.

## 5.2. Diocotron instability.
We now consider the diocotron instability [15] for an annular electron layer usually described by the guiding center model (1.7). This instability is well studied numerically as in [39, 31]. It may give rise to electron vortices, which is the analog of the Kelvin-Helmholtz fluid dynamics and may occur when charge neutrality is not locally maintained.

Here we want to investigate the development of such instability when we consider the Vlasov-Poisson system with an external magnetic field (1.3). It is expected that such instability holds true for large magnetic fields, whereas dissipative effects dominate when the magnetic field is not large enough to confine particles.

Therefore, we perform numerical simulation using a particle-in-cell method [5] for the Vlasov equation (1.3), where the particle trajectories are approximated thanks to the third order semi-implicit scheme (3.13)-(3.17). On the other hand, we compute an approximation of the guiding center model using a finite difference method developed in [39]. This reference solution will be used to compare our results obtained from the Vlasov-Poisson system with a large magnetic fields.

The initial density is given by

$$ \rho_0(\mathbf{x}) = \begin{cases} (1 + \alpha \cos(\ell\theta)) \exp(-4(\|\mathbf{x}\| - 6.5)^2), & \text{if } r^- \le \|\mathbf{x}\| \le r^+, \\ 0, & \text{otherwise,} \end{cases} $$

where $\alpha$ is a small parameter, $\theta = \operatorname{atan}(y/x)$, with $\mathbf{x} = (x, y) \in \mathbb{R}^2$. In the following tests, we take $\alpha = 0.01$, $r^- = 5$, $r^+ = 8$, $\ell = 7$.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Diocotron instability $\epsilon = 1$. Time evolution of (a) relative energy $\mathcal{E}(t) - \mathcal{E}(0)$ with respect to the initial one and (b) $\|\mathbf{E}(t)\|_{\infty}$. For $\epsilon \sim 1$ the instability does not occur.</figcaption></figure>

Furthermore, we will also consider the Vlasov-Poisson system (1.3) with an external magnetic field with the initial datum $f_0$

$$f_0(\mathbf{x}, \mathbf{v}) = \frac{\rho_0(\mathbf{x})}{2\pi} \exp\left(-\frac{\|\mathbf{v}\|^2}{2}\right), \quad (\mathbf{x}, \mathbf{v}) \in \mathbb{R}^4.$$

A particle method with a third order semi-implicit solver (3.13)-(3.17) will be applied for different values of $\epsilon = 1, 10^{-1}$ and $10^{-2}$.

For both systems a high order finite difference scheme in Cartesian coordinates will be applied [39] for the numerical approximation of the Poisson equation in a disc domain. We choose a grid with $100 \times 100$ points in the physical space and 400 particles per cell for the discretization of the velocity space.

We define the total energy $\mathcal{E}(t)$ as the sum of the potential energy $\mathcal{E}_{\text{pot}}(t)$ and the kinetic energy $\mathcal{E}_{\text{kin}}(t)$ with

$$\mathcal{E}_{\text{pot}}(t) = \frac{1}{2} \int_{\mathbb{R}^2} |\nabla \phi(t, \mathbf{x})|^2 d\mathbf{x}, \quad \text{and} \quad \mathcal{E}_{\text{kin}}(t) = \frac{1}{2} \int_{\mathbb{R}^4} f(t, \mathbf{x}, \mathbf{v}) |\mathbf{v}|^2 d\mathbf{x} d\mathbf{v}.$$

For the Vlasov-Poisson system (1.3), the total energy is exactly conserved with respect to time.

First, we consider the case where $\epsilon = 1$, that is, the particle trajectories do not coincide with the trajectory corresponding to the guiding center model (1.7). On the one hand, we present in Figures 3, the evolution of relative energy with respect to the initial data $t \mapsto \mathcal{E}(t) - \mathcal{E}(0)$ and the $L^{\infty}$ norm of the self-consistent electric field $t \mapsto \|\mathbf{E}(t)\|_{\infty}$. The discrete total energy is not exactly preserved but it oscillates around zero and these variations remain relatively small compared to the variations of the kinetic and potential energy which also oscillate with a frequency around $1/2\pi$. The same phenomenon can be observed on the time evolution of the $L^{\infty}$ norm of the electric field. On the other hand in Figure 4, we plot the time evolution of the charged density and observe that when the amplitude of the external magnetic field $\|\mathbf{B}\|$ is of order one, the plasma is not well confined and does not develop any instability. The density seems to oscillate around the steady state.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Diocotron instability $\varepsilon = 1$. Time evolution of the density $\rho$ for time $t = 36, t = 38, t = 40$ and $t = 42$ units. For $\varepsilon \sim 1$ the instability does not occur.</figcaption></figure>

Then, we take a small value $\varepsilon = 0.01$, and for this case we present a comparison between the finite difference approximation to the guiding center model (1.7) and a particle method for the Vlasov-Poisson system (1.3), using the third order semi-implicit scheme (3.13)-(3.17). The numerical results are presented in Figures 5 and 6.

Let us emphasize that for the guiding center model (1.7), that is in the limit $\varepsilon \to 0$, the potential energy $\mathcal{E}_{\text{pot}}$ is conserved with respect to time. Therefore, for $\varepsilon \ll 1$, it is expected that both the variations of $\mathcal{E}_{\text{pot}}$ and $\mathcal{E}_{\text{kin}}$ are small. In Figure 5, we can indeed observe that the variation of both quantities $\mathcal{E}_{\text{pot}}$ and $\mathcal{E}_{\text{kin}}$ are varying with a small amplitude of order $10^{-3}$. In that case, we do not see any oscillations since a large time step is used. Moreover, the second picture in Figure 5 represents the time evolution of $\|\mathbf{E}(\cdot)\|_{\infty}$ for the Vlasov-Poisson system (1.3) and the guiding center model (1.7). Both results agree well which illustrates the accuracy of the particle method in the limit $\varepsilon \ll 1$.
<figure><img src="image_6.png" /><figcaption>FIGURE 5. **Diocotron instability** $\epsilon = 10^{-2}$. Time evolution of the relative energy with respect to the initial one (a) for the Vlasov-Poisson system (1.3) (b) for the guiding center model (1.7).</figcaption></figure>

Finally in Figure 6, we present a comparison between the density $\rho$ obtained from the Vlasov-Poisson system discretized with the particle method and the one corresponding to the guiding center model discretized with a finite difference method. These figures show the development of the diocotron instability on the density $\rho$ for both models. Indeed in this regime ($\epsilon \ll 1$), it is expected that the density $\rho^\epsilon$ computed from the Vlasov-Poisson system obeys to the same evolution as the one of the guiding center model. Once again, the results agree very well and it is remarkable that the particle-in cell method does not suffer from too many fluctuations. The vortices are well described even for large time $t \sim 120$.

## 6. CONCLUSION AND PERSPECTIVE

In this paper we proposed a class of semi-implicit time discretization techniques for particle-in cell simulations. The main feature of this approach is to guarantee the accuracy and stability when the amplitude of the magnetic field becomes large and to get the correct long time behavior (guiding center approximation). We formally showed that the present schemes preserve the initial order of accuracy when $\epsilon \to 0$. Furthermore, we performed a complete analysis of the first order semi-implicit scheme when we consider a given and smooth electromagnetic field ($\mathbf{E}, \mathbf{B}_{\text{ext}}$).

The time discretization techniques proposed in this paper seem to be a very simple and efficient tool to filter fast oscillations and have nice stability and consistency properties in the limit $\epsilon \to 0$. However, a complete analysis of high order semi-implicit schemes is still missing. The main issue is to control the space trajectory $(\mathbf{x}^n)_n$ uniformly with respect to $\epsilon$ and as we have shown in Remark 4.2, the use of a semi-implicit scheme does not necessarily guarantee that the particle trajectories are under control. A complete analysis of high order schemes is currently under study.

On the other hand, the present techniques will be applied to more advanced problems as the three dimensional Vlasov-Poisson system when the magnetic field is non uniform and the particle trajectories become more complicated.
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Diocotron instability $\epsilon = 10^{-2}$. Time evolution of the density $\rho$ for time $t = 60$, $t = 90$, and $t = 120$ units for (a) the Vlasov-Poisson system (1.3) and (b) the guiding center model (1.7).</figcaption></figure>
ACKNOWLEDGEMENTS

FF was supported by the EUROfusion Consortium and has received funding from the Eu-
ratom research and training programme 2014-2018 under grant agreement No 633053. The
views and opinions expressed herein do not necessarily reflect those of the European Commis-
sion.

LMR was supported in part by the ANR project BoND (ANR-13-BS01-0009-01).

REFERENCES

[1] C. ALARD AND S. COLOMBI, A cloudy Vlasov solution. Monthly Notices of the Royal Astronomical Society, **359** (1) pp. 123163, (2005).

[2] W. B. BATESON AND D.W. HEWETT, Grid and Particle Hydrodynamics. Journal of Computational Physics, **144** pp. 358378, (1998).

[3] J. BEALE AND A. MAJDA, Vortex methods. II. Higher order accuracy in two and three dimensions. Mathematics of Computation **39** (159), pp. 29–52 (1982).

[4] B.A. DE DIOS, J.A. CARRILLO AND C.W. SHU, Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem, Mathematical Models and Methods in Applied Sciences, **22** (2012).

[5] C. K. BIRDSALL, A. B. LANGDON, Plasma Physics via Computer Simulation Institute of Physics Publishing, Bristol and Philadelphia., (1991).

[6] S. BOSCARINO, F. FILBET AND G. RUSSO, High order semi-Implicit schemes for time dependent partial differential equations, preprint (2014).

[7] M. CAMPOS PINTO, E. SONNENDRÜCKER, A. FRIEDMAN, D.P. GROTE, AND S.M. LUND, Noiseless Vlasov-Poisson simulations with linearly transformed particles. Journal of Computational Physics, **275** pp. 236256, (2014).

[8] J. CHENG, S. E. PARKER, Y. CHEN AND D. A. UZDENSKY, A second-order semi-implicit df method for hybrid simulation Journal of Computational Physics, **245** pp. 364375, (2013).

[9] A. COHEN AND B. PERTHAME, Optimal Approximations of Transport Equations by Particle and Pseudoparticle Methods. SIAM J. on Math. Anal. **32**(3), pp. 616–636 (2000).

[10] G.-H. COTTET AND P.-A. RAVIART, Particle methods for the one-dimensional Vlasov-Poisson equations, SIAM J. Numer. Anal., **21**, pp. 52-76 (1984).

[11] G.-H. COTTET AND P. KOUMOUTSAKOS, Vortex Methods: Theory and Practice. Cambridge University Press, Cambridge (2000)

[12] N. CROUSEILLES, E. FRÉNOD, S. HIRSTOAGA, A. MOUTON, Two-Scale Macro-Micro decomposition of the Vlasov equation with a strong magnetic field, Math. Models Methods Appl. Sci. **23** (2015)

[13] N. CROUSEILLES, M. MEHRENBERGER AND E. SONNENDRÜCKER, Conservative semi-Lagrangian schemes for Vlasov equations, Journal of Computational Physics, **229**, pp. 1927–1953 (2010).

[14] N. CROUSEILLES, M. LEMOU AND F. MÉHATS, Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations, Journal of Computational Physics, **248**, pp. 287–308 (2013).

[15] R. C. DAVIDSON, Physics of Nonneutral Plasmas, Imperial College Press, (2001).

[16] F. FILBET, Convergence of a Finite Volume Scheme for the One Dimensional Vlasov-Poisson System, *SIAM J. Numer. Analysis*, **39**, pp. 1146–1169 (2001).

[17] F. FILBET, E. SONNENDRÜCKER, P. BERTRAND, Conservative numerical schemes for the Vlasov equation, *Journal of Computational Physics*, **172**, pp. 166–187 (2001).

[18] F. FILBET AND E. SONNENDRÜCKER, Comparison of Eulerian Vlasov solvers, Computer Physics Communications, **150**, pp. 247–266 (2003).

[19] F. FILBET AND E. SONNENDRÜCKER, Modeling and numerical simulation of space charge dominated beams in the paraxial approximation, *Mathematical Models and Methods in the Applied Sciences*, **16**, pp. 763–791 (2006).

[20] E. FRÉNOD, S. HIRSTOAGA, M. LUTZ AND E. SONNENDRÜCKER, Long Time Behaviour of an Exponential Integrator for a Vlasov-Poisson System with Strong Magnetic Field. Commun. Comput. Phys. **18**, no. 2, 263296 (2015).

[21] E. FRÉNOD, F. SALVARINI, E. SONNENDRÜCKER, Long time simulation of a beam in a periodic focusing channel via a two-scale PIC-method, Math. Models Methods Appl. Sci. **10**, 175-197, (2009).

[22] E. FRÉNOD, E. SONNENDRÜCKER, Long time behavior of the Vlasov equation with strong external magnetic field, Math. Models Methods Appl. Sci. **19**, 539-553, (2000).
[23] K. GANGULY AND H. D. VICTORY, JR, On the convergence of particle methods for multidimensional Vlasov-Poisson systems, SIAM J. Numer. Anal., **26**, pp. 249–288 (1989).

[24] F. GOLSE AND L. SAINT RAYMOND, The Vlasov-Poisson system with strong magnetic field. J. Math. Pures. Appl., **78** pp. 791817, (1999).

[25] E. HAIRER, S.P. NORSETT AND G. WANNER, Solving ordinary differential equations. I. Non-stiff problems. Second edition. Springer Series in Computational Mathematics, 8. Springer-Verlag, Berlin, 1993. xvi+528 pp.

[26] R. E. HEATH, I. M. GAMBA, P.J. MORRISON AND C. MICHLER, A discontinuous Galerkin method for the Vlasov-Poisson system, *Journal of Computational Physics*, **231** (2012), pp. 1140–1174.

[27] D.W. HEWETT, Fragmentation, merging, and internal dynamics for PIC simulation with finite size particles. Journal of Computational Physics, **189** pp. 390426, (2003).

[28] S. JIN, Efficient asymptotic-preserving (AP) schemes for some multi scale kinetic equations, SIAM J. Sci. Comput. **21** pp. 441-454 (1999).

[29] A. KLAR, An asymptotic-induced scheme for non-stationary transport equations in the diffusive limit. SIAM J. Numer. Anal. **35**, no. 3, 10731094 (1998).

[30] P. KOUMOUTSAKOS, Inviscid Axisymmetrization of an Elliptical Vortex. Journal of Computational Physics **138**, 821–857 (1997).

[31] J. PÉTRI, Non-linear evolution of the diocotron instability in a pulsar electrosphere: 2D PIC simulations, *Astronomy & Astrophysics*, **503**, pp. 1–12 (2009).

[32] J. M. QIU AND C.-W. SHU, Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system, Journal of Computational Physics, **230**, pp. 8386–8409 (2011).

[33] P.-A. RAVIART, An analysis of particle methods. In: Numerical methods in fluid dynamics (Como, 1983), pp. 243–324. Lecture Notes in Mathematics, Berlin (1985).

[34] L. SAINT-RAYMOND, The gyro-kinetic approximation for the Vlasov-Poisson system. Math. Models Methods Appl. Sci. **10**, pp. 13051332 (2000).

[35] L. SAINT-RAYMOND, Control of large velocities in the two-dimensional gyro-kinetic approximation. J. Math. Pures Appl. **81**, no. 4, 379399 (2002).

[36] E. SONNENDRÜCKER, J. ROCHE, P. BERTRAND, A. GHIZZO, The semi-Lagrangian method for the numerical resolution of Vlasov equation, Journal of Computational Physics, **149**, pp. 201–220 (2013).

[37] S. WOLLMAN AND E. OZIZMIR, Numerical approximation of the one-dimensional Vlasov-Poisson system with periodic boundary conditions, SIAM J. Numer. Anal. **33**, pp. 1377–1409 (1996).

[38] S. WOLLMAN, On the approximation of the Vlasov-Poisson system by particle methods, SIAM J. Numer. Anal., **37**, pp. 1369–1398 (2000).

[39] C. YANG AND F. FILBET, Conservative and non-conservative methods based on Hermite weighted essentially non-oscillatory reconstruction for Vlasov equations. J. Comput. Phys. **279**, pp. 1836 (2014).

FRANCIS FILBET

UNIVERSITÉ DE TOULOUSE III & IUF

UMR5219, INSTITUT DE MATHÉMATIQUES DE TOULOUSE,

118, ROUTE DE NARBONNE

F-31062 TOULOUSE CEDEX, FRANCE

E-MAIL: francis.filbet@math.univ-toulouse.fr
LUIS MIGUEL RODRIGUES

UNIVERSITÉ DE RENNES 1,
UMR6625, IRMAR,
263 AVENUE DU GENERAL LECLERC,
F-35042 RENNES CEDEX, FRANCE

E-MAIL: luis-miguel.rodrigues@univ-rennes1.fr