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

HAL
open science

# A modified Crank-Nicolson scheme for the Vlasov-Poisson system with a strong external magnetic field

Francis Filbet, L Miguel Rodrigues, Kim Han Trinh

► To cite this version:

Francis Filbet, L Miguel Rodrigues, Kim Han Trinh. A modified Crank-Nicolson scheme for the Vlasov-Poisson system with a strong external magnetic field. 2025. hal-05138815v3

HAL Id: hal-05138815

https://hal.science/hal-05138815v3

Preprint submitted on 2 Dec 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.

Public Domain
# A MODIFIED CRANK-NICOLSON SCHEME FOR THE VLASOV-POISSON SYSTEM WITH A STRONG EXTERNAL MAGNETIC FIELD

FRANCIS FILBET, L. MIGUEL RODRIGUES, AND KIM HAN TRINH

**ABSTRACT.** We propose and study a Particle-In-Cell (PIC) method utilizing Crank-Nicolson time discretization for the Vlasov-Poisson system with a strong, inhomogeneous external magnetic field with fixed direction. Our focus is on particle dynamics in the plane orthogonal to the magnetic field. In this regime, traditional explicit schemes are constrained by stability conditions linked to the small Larmor radius and plasma frequency [23]. To avoid this limitation, our approach is based on numerical schemes [11, 12, 14], providing a consistent PIC discretization of the guiding-center system taking into account variations of the magnetic field. We carry out some theoretical proofs and perform several numerical experiments to validate the method demonstrating its robustness and accuracy.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>1.1. Formal asymptotic behavior for a given electromagnetic field</td><td>2</td></tr><tr><td>1.2. Formal asymptotic limit of the Vlasov-Poisson system</td><td>3</td></tr><tr><td>1.3. Particle methods for the Vlasov-Poisson system</td><td>4</td></tr><tr><td>2. Review of Crank-Nicolson schemes</td><td>5</td></tr><tr><td>2.1. The Brackbill-Forslund-Vu scheme</td><td>8</td></tr><tr><td>2.2. The Ricketson-Chacón scheme</td><td>10</td></tr><tr><td>3. A modified Crank-Nicolson scheme with an additional variable</td><td>12</td></tr><tr><td>4. Numerical simulations</td><td>13</td></tr><tr><td>4.1. One single particle motion</td><td>14</td></tr><tr><td>4.2. Vlasov-Poisson system</td><td>18</td></tr><tr><td>5. Conclusion</td><td>23</td></tr><tr><td>Appendix A. Formal asymptotic behavior for a given electromagnetic field</td><td>25</td></tr><tr><td>References</td><td>28</td></tr></table>

## 1. INTRODUCTION

This paper focuses on plasma confinement in the presence of a strong, spatially varying external magnetic field, where charged particles evolve under the influence of both electrostatic forces and intense magnetic confinement. Such configurations are characteristic of tokamak plasmas [1, 21], where the magnetic field plays a crucial role in containing particles within the core of the device. Kinetic models, which provide a mesoscopic description of charged particle dynamics, are highly accurate and essential tools for investigating the behavior of thermonuclear fusion plasmas.

We assume that collective effects dominate, and the plasma is modeled entirely through **kinetic** equations. The primary unknown is the particle number density $f \equiv f(t, \mathbf{x}, \mathbf{v})$, which depends on time $t \ge 0$, position $\mathbf{x} \in \Omega \subset \mathbb{R}^d$, and velocity $\mathbf{v} \in \mathbb{R}^d$, with $d \ge 2$. 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 $F(t, \mathbf{x}, \mathbf{v})$ is coupled with the distribution function $f$ giving a nonlinear system.

<footnote n="2010">2010 Mathematics Subject Classification. Primary: 65M75 Secondary: 82D10 76X05 35Q83.</footnote>

<footnote>Key words and phrases. Vlasov-Poisson systems; Strong magnetic field; Particle methods.</footnote>

<footnote>Research of L.M.R. was partially supported by the ANR Project HEAD ANR-24-CE40-3260 and the Institut Universitaire de France.</footnote>

<footnote>Research of K.H.T. was partially supported by the Institut Universitaire de France.</footnote>
Here, we consider the two-dimensional case where the magnetic field acts in the vertical direction and only depends on $x = (x_1, x_2) \in \mathbb{R}^2$, that is,

$$ \mathbf{B}(\mathbf{x}) = \frac{1}{\varepsilon} \begin{pmatrix} 0 \\ 0 \\ b(\mathbf{x}) \end{pmatrix}, $$

where the function $b$ describes the variations of the amplitude with $b \in W^{1,\infty}(\mathbb{R}^2)$ and

$$ (1.2) \qquad b(\mathbf{x}) \geq b_0 > 0. $$

The number $\varepsilon > 0$ is a small parameter related to the ratio between the reciprocal Larmor frequency and the advection time scale (see [10, 19, 20] and the references therein for more details on scalings).

We will focus on the long-time behavior of positive ions in the orthogonal plane to the external magnetic field. Therefore, the distribution function $f_\varepsilon$ is a solution to the Vlasov equation coupled with the Poisson equation for the electrical potential $\phi_\varepsilon$ generated by the motion of these charged particles, that is,

$$ (1.3) \qquad \begin{cases} \varepsilon \frac{\partial f_\varepsilon}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f_\varepsilon + \left( \mathbf{E}_\varepsilon(t, \mathbf{x}) - b(\mathbf{x}) \frac{\mathbf{v}^\perp}{\varepsilon} \right) \cdot \nabla_{\mathbf{v}} f_\varepsilon = 0, \\ \mathbf{E}_\varepsilon = -\nabla_{\mathbf{x}} \phi_\varepsilon, \quad -\Delta_{\mathbf{x}} \phi_\varepsilon = \rho_\varepsilon, \end{cases} $$

where $\mathbf{v}^\perp = (-v_2, v_1) \in \mathbb{R}^2$ and the density $\rho_\varepsilon$ is given by

$$ \rho_\varepsilon(t, \mathbf{x}) := \int_{\mathbb{R}^2} f_\varepsilon(t, \mathbf{x}, \mathbf{v}) \, d\mathbf{v}. $$

Here we aim to construct numerical approximations for the Vlasov-Poisson system (1.3) using particle methods (see [2]), which involve in approximating the distribution function by a finite number of macro-particles. The trajectories of these particles are determined from the characteristic curves associated to the Vlasov equation

$$ (1.4) \qquad \left\{ \begin{aligned} \varepsilon \frac{d\mathbf{x}_\varepsilon}{dt} &= \mathbf{v}_\varepsilon, \\ \varepsilon \frac{d\mathbf{v}_\varepsilon}{dt} &= \mathbf{E}_\varepsilon(t, \mathbf{x}_\varepsilon) - b(\mathbf{x}_\varepsilon) \frac{\mathbf{v}_\varepsilon^\perp}{\varepsilon}, \\ \mathbf{x}_\varepsilon(0) &= \mathbf{x}_\varepsilon^0, \quad \mathbf{v}_\varepsilon(0) = \mathbf{v}_\varepsilon^0, \end{aligned} \right. $$

then we use the conservation of $f_\varepsilon$ along the characteristic curves, that is,

$$ f_\varepsilon(t, \mathbf{x}_\varepsilon(t), \mathbf{v}_\varepsilon(t)) = f_\varepsilon(t^0, \mathbf{x}_\varepsilon^0, \mathbf{v}_\varepsilon^0). $$

In particular, we will focus on the construction of numerical schemes for the ODE system (1.4), where the time step $\Delta t$ is arbitrary and free from any stability constraint. Following the work of Filbet and Rodrigues [12–14], the ODE system can be decomposed into fast dynamics, driven by the fast variable $\mathbf{v}_\varepsilon$, and slow dynamics, governed by the variables $(\mathbf{x}_\varepsilon, e_\varepsilon)$, where $e_\varepsilon = \frac{1}{2}|\mathbf{v}_\varepsilon|^2$. This decomposition allows the design of a class of numerical schemes that precisely capture slow-scale variables, while faster scales are correctly filtered. More precisely, when the intensity of the magnetic field is sufficiently large, i.e. when $\varepsilon \ll 1$, the scheme provides a consistent approximation to the asymptotic model [12].

**1.1. Formal asymptotic behavior for a given electromagnetic field.** Before describing a numerical scheme for the nonlinear Vlasov-Poisson system (1.3), we first briefly expound on what may be expected from the continuous model with a given electric field in the limit $\varepsilon \to 0$. For this purpose, we consider a function $\phi \in W^{3,\infty}$ such that, for all $\mathbf{x} \in \mathbb{R}^2$, $\mathbf{E}(\mathbf{x}) = -\nabla_{\mathbf{x}}\phi(\mathbf{x})$ and observe that the system (1.4) has an Hamiltonian structure associated with the total energy $\mathcal{E}_\varepsilon(t)$,

$$ (1.5) \qquad \mathcal{E}_\varepsilon(t) = \frac{\|\mathbf{v}_\varepsilon(t)\|^2}{2} + \phi(\mathbf{x}_\varepsilon(t)), \quad t \ge 0, $$

which is an invariant of the system. Therefore, to study the limit $\varepsilon \to 0$, we first define the kinetic energy as a slow scale variable

$$ e_\varepsilon(t) := \frac{1}{2} \|\mathbf{v}_\varepsilon(t)\|^2, $$
leading to the study of the augmented system

$$ (1.6) \qquad \begin{cases} \frac{d\mathbf{x}_\varepsilon}{dt} = \frac{\mathbf{v}_\varepsilon}{\varepsilon}, \\ \frac{de_\varepsilon}{dt} = \frac{1}{\varepsilon} \mathbf{E}(\mathbf{x}_\varepsilon) \cdot \mathbf{v}_\varepsilon, \end{cases} $$

still coupled with the equation on $\mathbf{v}_\varepsilon$

$$ (1.7) \qquad \varepsilon \frac{d\mathbf{v}_\varepsilon}{dt} = \mathbf{E}(\mathbf{x}_\varepsilon) - b(\mathbf{x}_\varepsilon) \frac{\mathbf{v}_\varepsilon^\perp}{\varepsilon}, $$

which describes the fastest scale. Of course the second equation of (1.6) is a consequence of (1.7) but it retains only its slower part.

Following [13], one may prove that $\mathbf{x}_\varepsilon(t) \to \mathbf{y}(t)$ and $e_\varepsilon(t) \to g(t)$, as $\varepsilon \to 0$ where $(\mathbf{y}, g)$ is solution to the so-called guiding center system,

$$ (1.8) \qquad \begin{cases} \frac{d\mathbf{y}}{dt} = -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}) + g \frac{\nabla_y^\perp b}{b^2}(\mathbf{y}), \\ \frac{dg}{dt} = g \mathbf{E} \cdot \frac{\nabla_y^\perp b}{b^2}(\mathbf{y}). \end{cases} $$

For the convenience of the reader we provide in Appendix A the main formal computations leading to (1.8).

Furthermore, we may identify the limit for the total energy $\mathcal{E}_\varepsilon(t)$, as $\varepsilon \to 0$, this quantity converging to

$$ \mathcal{E}_{gc}(t) := g(t) + \phi(\mathbf{y}(t)), $$

which is indeed an invariant of the guiding center model (1.8),

$$ (1.9) \qquad \frac{d\mathcal{E}_{gc}}{dt} = \frac{dg}{dt} + \nabla_y \phi(\mathbf{y}) \cdot \frac{d\mathbf{y}}{dt} = g \mathbf{E} \cdot \frac{\nabla_y^\perp b}{b^2}(\mathbf{y}) - g \mathbf{E} \cdot \frac{\nabla_y^\perp b}{b^2}(\mathbf{y}) = 0. $$

Furthermore, we may define the magnetic moment as

$$ \mu_{gc}(t) = \frac{g}{b(\mathbf{y})}, $$

which is an invariant for the guiding center system (1.8) without counterpart for the original (1.4), thus called an adiabatic invariant for (1.8). Indeed, we have

$$ (1.10) \qquad \frac{d\mu_{gc}}{dt} = \frac{d}{dt} \left( \frac{g}{b(\mathbf{y})} \right) = g \mathbf{E} \cdot \frac{\nabla_y^\perp b}{b^3}(\mathbf{y}) - g \frac{\nabla_y b}{b^2}(\mathbf{y}) \cdot \left( -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}) + g \frac{\nabla_y^\perp b}{b^2}(\mathbf{y}) \right) = 0. $$

Reproducing these properties at the discrete level is a target when designing a scheme for (1.4) preserving asymptotics when $\varepsilon$ tends to zero.

1.2. **Formal asymptotic limit of the Vlasov-Poisson system.** We come back to the Vlasov-Poisson system (1.3). Here one cannot anymore remain completely at the characteristic level (1.4). Moreover whereas arguments of the previous subsection could be turned into sound analytic arguments (by slight variations on [13]), to the best of our knowledge the present situation does not fall directly into the range of the actually available analysis of gyro-kinetic limits. We refer the reader to the introductions of [13, 25] and references therein for a representative sample of such analytic techniques.

Nevertheless the known results and the previous subsection strongly suggests for $(f^\varepsilon, \mathbf{E}^\varepsilon)$ solving the Vlasov-Poisson system (1.3) that in the limit $\varepsilon \to 0$, the electric field $\mathbf{E}^\varepsilon$ and the following velocity-averaged version of $\bar{F}^\varepsilon$

$$ \bar{F}^\varepsilon : (t, \mathbf{x}, e) \mapsto \frac{1}{2\pi} \int_0^{2\pi} f^\varepsilon(t, \mathbf{x}, \sqrt{2e}(\cos(\theta), \sin(\theta))) d\theta $$
converge to some $\mathbf{E} : (t, \mathbf{y}) \mapsto \mathbf{E}(t, \mathbf{y})$ and some$^1$ $f : (t, \mathbf{y}, g) \mapsto f(t, \mathbf{y}, g)$ solving the following system consisting in a kinetic equation supplemented with a Poisson equation,

$$ (1.11) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial f}{\partial t} + \mathbf{U} \cdot \nabla_{\mathbf{y}} f + u_g \frac{\partial f}{\partial g} = 0, \\ \\ -\Delta_{\mathbf{y}} \phi = \rho, \quad \rho = 2\pi \int_{\mathbb{R}^+} f \, dg, \end{array} \right. $$

where the velocity field is given by

$$ \mathbf{U}(t, \mathbf{y}, g) = \mathbf{F}(t, \mathbf{y}) + g \frac{\nabla_{\mathbf{y}} b}{b^2}(t, \mathbf{y}), \quad u_g = -\mathrm{div}_{\mathbf{y}}(\mathbf{F})(t, \mathbf{y}) g, $$

with $\mathbf{E} = -\nabla_{\mathbf{y}}\phi$, $\mathbf{F} = -\mathbf{E}^{\perp}/b$. We remind the reader that $\mathbf{U}$ contains two classical components of the guiding center velocity, the $\mathbf{E} \times \mathbf{B}$ drift and the grad $\mathbf{B}$ drift.

### 1.3. Particle methods for the Vlasov-Poisson system.
To make the most of the previous discussions in order to discretize the Vlasov equation (1.3), particle methods are particularly well suited since they directly involve an approximation of the characteristic curves (1.4). Here, we will consider the Particle-In-Cell (PIC) method, in which trajectories are computed via the characteristic curves (1.4), while the self-consistent electric field is calculated using the Poisson equation on a grid of the physical space. We refer the reader to [2,9] or [11] for a brief review of particle methods.

To keep the notation as concise as possible, we temporarily omit the dependence of solutions on $\varepsilon$. The starting point is the approximation of the solution $f$, which solves (1.3), by a finite sum of smoothed functions — viewed as macro particles. More explicitly, in dimension $d$, one computes

$$ f_N(t, \mathbf{x}, \mathbf{v}) = \sum_{1 \le k \le N} \omega_k \varphi_\alpha (\mathbf{x} - \mathbf{x}_k(t)) \otimes \varphi_\alpha (\mathbf{v} - \mathbf{v}_k(t)), $$

where $\varphi_\alpha = \alpha^{-d}\varphi(\cdot/\alpha)$ is a particle shape function with radius proportional to $\alpha$ — usually seen as an approximation of the Dirac measure $\delta_0$ — obtained by rescaling a fixed compactly supported mollifier $\varphi$ whereas the set $(\mathbf{x}_k, \mathbf{v}_k)_{1 \le k \le N}$ represents the position in phase space of $N$ macro-particles evolving along characteristic curves (1.4) from the initial data $(\mathbf{x}_k^0, \mathbf{v}_k^0)$, $1 \le k \le N$. More explicitly, $(\mathbf{x}_k, \mathbf{v}_k)_{1 \le k \le N}$ is solution to

$$ \left\{ \begin{array}{l} \varepsilon \frac{d\mathbf{x}_k}{dt} = \mathbf{v}_k, \\ \\ \varepsilon \frac{d\mathbf{v}_k}{dt} = \mathbf{E}(t, \mathbf{x}_k) - b(t, \mathbf{x}_k) \frac{\mathbf{v}_k^{\perp}}{\varepsilon}, \\ \\ \mathbf{x}_k(0) = \mathbf{x}_k^0, \quad \mathbf{v}_k(0) = \mathbf{v}_k^0, \end{array} \right. $$

where the electric field $\mathbf{E}$ is computed by discretizing the Poisson equation on a mesh of the physical space.

Here, we deliberately choose to use a classical PIC method with $P_1$ reconstruction for the density and electric field in order to focus on resolving particle trajectories with a time step large compared to the $\varepsilon$ scale parameter. Our approach can be easily extended to much more sophisticated PIC methods. For example, forward-backward Lagrangian methods [6] reduce fluctuations in the density reconstruction step. In recent years, other advanced methods have been developed to improve the stability properties of the Particle-In-Cell (PIC) method$^2$ in the presence of a large, inhomogeneous external magnetic field. Among these, the earliest schemes were introduced by Boris [2,3] for relativistic plasma simulation. It is a second-order explicit method, often referred to as an explicit PIC method, employing a time-centered electromagnetic field and an averaged phase-space representation $(\mathbf{x}, \mathbf{v})$ for (1.4). Later, this scheme was applied by Parker and Birdsall [22] to address the high magnetic field regime, aiming to accurately capture drift motions of particles in three dimensions. However, these standard explicit PIC approaches still suffer from temporal numerical stability constraints imposed by the Courant-Friedrichs-Lewy (CFL) condition [2]. As a result, despite their simplicity and computational efficiency, these schemes are significantly constrained in high-field regimes.

<footnote>$^1$We use distinct notation of variables for limiting functions to be consistent with asymptotic analysis at the characteristic level. This is of course completely immaterial.</footnote>

<footnote>$^2$For a discussion of some other classes of methods, we refer to the introduction of [14] and some references therein.</footnote>
To overcome this lack of stability, several implicit PIC schemes have been developed to solve (1.4) and to capture grad **B** drifts in strongly magnetized plasma. We refer to Brackbill, Forslund, and Vu [5, 26], who introduce an effective force into the velocity equation such that the scheme remains consistent with (1.4) for small time steps. The scheme is formulated in a fully implicit manner as a modified version of the Crank-Nicolson scheme. Alternatively, the magnetized implicit (MI) scheme proposed by Genoni, Clark, and Welch [18] employs a two-step predictor-corrector approximation. However, these schemes overlook the role of kinetic energy, which significantly contributes to particle motion when $\epsilon \to 0$, as shown in (1.8). Consequently, these schemes fail to capture the correct regime when, for a fixed time step, $\epsilon \to 0$.

More recently, Ricketson, Chacón and Chen [7, 23] built upon the Crank-Nicolson scheme with an additional effective force designed to achieve two objectives. First, they addressed the challenges arising in the regime $\epsilon \ll 1$, by capturing grad **B** drifts. Second, the additional force is designed to conserve energy for all $\epsilon > 0$. However, this method still requires to adapt the time step when $\epsilon$ becomes small. In parallel, building on [4], Filbet and Rodrigues proposed a class of semi-implicit methods (IMEX) where the position is updated explicitly whereas velocity is treated implicitly [11, 12, 14]. These schemes are developed to solve the augmented system (1.6), which introduces additional variable to separate slow scale and fast scale dynamics. As a result, in the regime $\epsilon \to 0$ and for a fixed time step, IMEX schemes accurately describe the dynamics of both position and kinetic energy. This ensures a consistent approximation to the guiding center system (1.6).

Another strategy has been developed to accurately capture the dynamics of (1.4) by following fast oscillations. This approach works well when the magnetic field is constant or when it varies slowly. For instance, the two-scaled formulation method, proposed in [8], employs two time variables to split fast and slow scales. Additionally, a class of Lie-Trotter type splitting schemes coupled with exponential integrators has been developed by Wang and Zhao [27] to provide a first order approximation with respect to $\epsilon$. These schemes are very successful for uniform or slowly varying magnetic field, but require a deeply understanding of the fast oscillations.

In this article, we propose to delve deeper and extend the strategy already proposed by Filbet and Rodrigues [11, 12] for Crank-Nicolson-type schemes. On the one hand, these schemes are widely recognized in the computational physics community and are valued for their effective energy preservation [5, 23, 26]. On the other hand, it is important to note that the schemes proposed by Filbet and Rodrigues, relying on IMEX methods that are more dissipative, can sometimes compromise their accuracy for intermediate values of the parameter $\epsilon$. The present work aims to maximize their efficiency and robustness by applying Crank-Nicolson schemes to the augmented system (1.6) to separate slow and fast scale dynamics. Moreover, this numerical scheme is implemented within the PIC framework for long-term simulations of the Vlasov-Poisson system (1.3).

The rest of the paper is organized as follows. In Section 2, we recall and analyze several numerical schemes based on the Crank-Nicolson scheme, including those developed by Brackbill, Forslund, and Vu [5, 26] as well as Ricketson and Chacón [23]. These schemes will be analyzed in the asymptotic limit $\epsilon \to 0$ with a fixed time step to clarify the importance of decomposing the solution into fast and slow variables, as in [11, 12]. In Section 3, we propose a new numerical scheme built upon the Crank-Nicolson scheme, following the strategy of Filbet and Rodrigues in [12], and we examine its accuracy in the regime $\epsilon \to 0$. Finally, in Section 4, we present numerical experiments for the new scheme, both for the computation of single-particle motion and as a particle pusher within the PIC framework for the Vlasov-Poisson system.

**Acknowledgement.** KHT expresses his appreciation of the hospitality of IMT, Université Toulouse III, during the preparation of the present contribution. FF and LMR are grateful to Luis Chacón for stimulating discussions that have motivated the present piece of work.

## 2. REVIEW OF CRANK-NICOLSON SCHEMES

In this section, we aim to discuss the Crank-Nicolson method, applied in the framework of Particle-In-Cell methods for the Vlasov-Poisson system. In [24], the authors show that the Crank-Nicolson scheme is second-order accurate, unconditionally stable, and energy-conserving for quadratic potentials when considering the system for a single particle motion (1.4). Later, this scheme has been studied for the Vlasov-Poisson system [7, 23, 26].

Here, we will review different modified Crank-Nicolson schemes proposed in the literature and discuss their conservation properties and asymptotic behavior when $\epsilon$ approaches zero, that is, when the external
magnetic field becomes large. Our aim is to investigate the consistency of the numerical approximation with the guiding center model (1.8) at the discrete level in the limit as $\epsilon \to 0$.

Let us start with the original Crank-Nicolson scheme and consider a time step $\Delta t > 0$ and $t^n = n \Delta t$, for $n \in \mathbb{N}$, we define $(\mathbf{x}_\epsilon^n, \mathbf{v}_\epsilon^n)$ as an approximation of the solution $(\mathbf{x}_\epsilon, \mathbf{v}_\epsilon)$ to (1.4). Applying the Crank-Nicolson scheme, the sequence $(\mathbf{x}_\epsilon^n, \mathbf{v}_\epsilon^n)_{n \in \mathbb{N}}$ is given by

$$ (2.1) \qquad \left\{ \begin{array}{l} \epsilon \frac{\mathbf{x}_\epsilon^{n+1} - \mathbf{x}_\epsilon^n}{\Delta t} = \mathbf{v}_\epsilon^{n+1/2}, \\ \epsilon \frac{\mathbf{v}_\epsilon^{n+1} - \mathbf{v}_\epsilon^n}{\Delta t} = \mathbf{E}(\mathbf{x}_\epsilon^{n+1/2}) - b(\mathbf{x}_\epsilon^{n+1/2}) \frac{(\mathbf{v}_\epsilon^{n+1/2})^\perp}{\epsilon}, \end{array} \right. $$

where

$$ \mathbf{v}_\epsilon^{n+1/2} = \frac{\mathbf{v}_\epsilon^{n+1} + \mathbf{v}_\epsilon^n}{2} \quad \text{and} \quad \mathbf{x}_\epsilon^{n+1/2} = \frac{\mathbf{x}_\epsilon^{n+1} + \mathbf{x}_\epsilon^n}{2}. $$

First, it is worth mentioning that this scheme provides a good approximation of the energy for a wide range of values of $\epsilon$. To be precise, assume that the electric field derives from a given smooth potential $\phi$, hence we have $\mathbf{E} = -\nabla_x \phi$ and the discrete kinetic energy is defined as

$$ e_{\epsilon}^{n} = \frac{1}{2} \|\mathbf{v}_{\epsilon}^{n}\|^{2}. $$

The total discrete energy is then

$$ (2.2) \qquad \mathcal{E}_{\epsilon}^{n} = e_{\epsilon}^{n} + \phi(\mathbf{x}_{\epsilon}^{n}), \quad n \ge 0. $$

From (2.1), it follows that the variation of the total energy is given by

$$ (2.3) \qquad \begin{aligned} \frac{\mathcal{E}_{\epsilon}^{n+1} - \mathcal{E}_{\epsilon}^{n}}{\Delta t} &= \frac{e_{\epsilon}^{n+1} - e_{\epsilon}^{n}}{\Delta t} + \frac{\phi(\mathbf{x}_{\epsilon}^{n+1}) - \phi(\mathbf{x}_{\epsilon}^{n})}{\Delta t}, \\ &= -\nabla_{\mathbf{x}}\phi(\mathbf{x}_{\epsilon}^{n+1/2}) \cdot \frac{\mathbf{x}_{\epsilon}^{n+1} - \mathbf{x}_{\epsilon}^{n}}{\Delta t} + \frac{\phi(\mathbf{x}_{\epsilon}^{n+1}) - \phi(\mathbf{x}_{\epsilon}^{n})}{\Delta t}. \end{aligned} $$

where $\mathbf{E}(\mathbf{x}_\epsilon^{n+1/2}) = -\nabla_\mathbf{x}\phi(\mathbf{x}_\epsilon^{n+1/2})$. Therefore, as is well-known, this scheme conserves the discrete energy only for quadratic potentials and for more general potential $\phi \in W^{3,\infty}$, a Taylor expansion yields

$$ (2.4) \qquad \phi(\mathbf{x}_{\epsilon}^{n+1}) - \phi(\mathbf{x}_{\epsilon}^{n}) = \nabla_{\mathbf{x}}\phi(\mathbf{x}_{\epsilon}^{n+1/2}) \cdot (\mathbf{x}_{\epsilon}^{n+1} - \mathbf{x}_{\epsilon}^{n}) + \Delta t^3 \mathcal{O}\left(\left\|\frac{\mathbf{x}_{\epsilon}^{n+1} - \mathbf{x}_{\epsilon}^{n}}{\Delta t}\right\|^3\right). $$

thus

$$ (2.5) \qquad \frac{\mathcal{E}_{\epsilon}^{n+1} - \mathcal{E}_{\epsilon}^{n}}{\Delta t} = \Delta t^2 \mathcal{O}\left(\left\|\frac{\mathbf{x}_{\epsilon}^{n+1} - \mathbf{x}_{\epsilon}^{n}}{\Delta t}\right\|^3\right). $$

In other words, under our above assumptions and the further reasonable assumption that $\mathbf{x}_\epsilon^n$ is bounded, the variations of the total discrete energy is of order $\Delta t^2$, which endows this scheme with a form of stability for large time simulations and for all $\epsilon > 0$.

Now, concerning the asymptotic limit of the scheme (2.1) as $\epsilon$ tends to zero, we have the following result.

**Proposition 2.1** (Asymptotic behavior $\epsilon \to 0$ with a fixed $\Delta t$). Let $\phi \in W^{3,\infty}(\mathbb{R}^2)$, choose a sufficiently small fixed time step $\Delta t$ and a final time $T > 0$. We set $N_T = [T/\Delta t]$. Assume that the Crank-Nicolson scheme (2.1) defines a numerical approximation $(\mathbf{x}_\epsilon^n, \mathbf{v}_\epsilon^n)_{0 \le n \le N_T}$ satisfying

(i) for all $1 \le n \le N_T$, $\mathbf{x}_\epsilon^n$ is uniformly bounded with respect to $\epsilon > 0$;

(ii) in the limit $\epsilon \to 0$, $(\mathbf{x}_\epsilon^0, \frac{1}{2}\|\mathbf{v}_\epsilon^0\|^2)$ converges to some $(\mathbf{y}^0, g^0)$.

Then, we have

* for all $1 \le n \le N_T$, $(\mathbf{x}_\epsilon^n, e_\epsilon^n)$ converges to $(\mathbf{y}^n, g^n)$, as $\epsilon \to 0$ with $e_\epsilon^n = \frac{1}{2}\|\mathbf{v}_\epsilon^n\|^2$ and the limit $(\mathbf{y}^n, g^n)_{1 \le n \le N_T}$ solves

$$ (2.6) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}^{n+1/2}), \\ \\ \displaystyle \frac{g^{n+1} - g^n}{\Delta t} = 0; \end{array} \right. $$
• for all $1 \le n \le N_T$, the total energy $\mathcal{E}_\varepsilon^n = e_\varepsilon^n + \phi(\mathbf{x}_\varepsilon^n)$ converges to $\mathcal{E}_{gc}^n := g^n + \phi(\mathbf{y}^n)$ as $\varepsilon \to 0$, which satisfies

$$ (2.7) \qquad \frac{\mathcal{E}_{gc}^{n+1} - \mathcal{E}_{gc}^n}{\Delta t} = \mathcal{O}(\Delta t^2); $$

• defining the discrete magnetic moment $\mu_\varepsilon^n$ as

$$ \mu_\varepsilon^n = \frac{e_\varepsilon^n}{b(\mathbf{x}_\varepsilon^n)}, $$

$(\mu_\varepsilon^n)_{\varepsilon>0}$ converges to $\mu_{gc}^n := \frac{g^n}{b(\mathbf{y}^n)}$ as $\varepsilon \to 0$ such that

$$ (2.8) \qquad \frac{\mu_{gc}^{n+1} - \mu_{gc}^n}{\Delta t} = -g^0 \frac{\mathbf{E} \cdot \nabla_{\mathbf{y}}^{\perp} b}{b^3} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2), $$

where $\mathbf{y}^{n+1/2}$ is defined as $\mathbf{y}^{n+1/2} = (\mathbf{y}^{n+1} + \mathbf{y}^n)/2$.

*Proof.* From our first assumption and the first equation of (2.1), we derive that each $(\varepsilon^{-1}\mathbf{v}_\varepsilon^{n+1/2})_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon$. By taking the limit $\varepsilon \to 0$ in the triangle inequality $\|\mathbf{v}_\varepsilon^{n+1}\| - \|\mathbf{v}_\varepsilon^n\| \le 2\|\mathbf{v}_\varepsilon^{n+1/2}\|$, one then deduces the convergence of $e_\varepsilon^n$ and the second equation of (2.6), thus also a uniform bound on $(\mathbf{v}_\varepsilon^n)_{\varepsilon>0}$.

Now, let us extract a subsequence still abusively labeled by $\varepsilon$ such that $\mathbf{x}_\varepsilon^n$ converges to some $\mathbf{y}^n$ as $\varepsilon$ goes to zero. By using the derived boundedness one may then take first a limit in the second equation of (2.1) to obtain

$$ \lim_{\varepsilon \to 0} \frac{\mathbf{v}_\varepsilon^{n+1/2}}{\varepsilon} = -\frac{1}{b(\mathbf{y}^{n+1/2})} \mathbf{E}^{\perp}(\mathbf{y}^{n+1/2}). $$

Then, take a limit in the first equation of (2.1) to conclude the derivation of (2.6). The latter uniquely characterizes the limit of the subsequence, thereby implying full convergence. At this stage, we use the smallness of $\Delta t$ (independent of $T$ and $\varepsilon$) to guarantee that the implicit scheme (2.6) is indeed solvable.

Let us now turn to the evolution of the total energy $\mathcal{E}_\varepsilon^n$ and the magnetic moment $\mu_\varepsilon^n$. For any $0 \le n \le N_T - 1$, the convergence of $(\mathbf{x}_\varepsilon^n, e_\varepsilon^n)$ to $(\mathbf{y}^n, g^n)$ as $\varepsilon$ goes to zero implies the convergence of $(\mathcal{E}_\varepsilon^n, \mu_\varepsilon^n)$ to $(\mathcal{E}_{gc}^n, \mu_{gc}^n)$. The total energy of the limiting system $\mathcal{E}_{gc}^n$ satisfies

$$ \frac{\mathcal{E}_{gc}^{n+1} - \mathcal{E}_{gc}^n}{\Delta t} = \frac{g^{n+1} - g^n}{\Delta t} + \frac{\phi(\mathbf{y}^{n+1}) - \phi(\mathbf{y}^n)}{\Delta t} = -\nabla_{\mathbf{y}}\phi(\mathbf{y}^{n+1/2}) \cdot \frac{(\mathbf{y}^{n+1} - \mathbf{y}^n)}{\Delta t} + \mathcal{O}(\Delta t^2) = \mathcal{O}(\Delta t^2), $$

as deduced from a Taylor expansion and the insertion of $\mathbf{E}(\mathbf{y}^{n+1/2}) = -\nabla_{\mathbf{y}}\phi(\mathbf{y}^{n+1/2})$ in the first equation of (2.6). Similarly, the evolution of the discrete magnetic moment $\mu_{gc}^n$ obeys for all $0 \le n \le N_T - 1$,

$$ \begin{aligned} \frac{\mu_{gc}^{n+1} - \mu_{gc}^n}{\Delta t} &= \frac{1}{\Delta t} \left( \frac{g^{n+1}}{b(\mathbf{y}^{n+1})} - \frac{g^n}{b(\mathbf{y}^n)} \right) \\ &= -g^0 \frac{\nabla_{\mathbf{y}} b}{b^2} (\mathbf{y}^{n+1/2}) \cdot \frac{(\mathbf{y}^{n+1} - \mathbf{y}^n)}{\Delta t} + \mathcal{O}(\Delta t^2) \\ &= -g^0 \frac{\mathbf{E} \cdot \nabla_{\mathbf{y}}^{\perp} b}{b^3} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2). \end{aligned} $$

□

It is worth mentioning that Proposition 2.1, clearly indicates that as $\varepsilon$ goes to zero, the discrete guiding center system (2.6) obtained by passing to the limit in the *Crank-Nicolson scheme* is not consistent with the continuous system (1.8). Indeed, it fails to capture the correct drift $\nabla_{\mathbf{y}}^{\perp} b/b^2$ for both position $\mathbf{y}$ and kinetic energy $g$. Similarly, the evolution of the magnetic moment $(\mu_{gc}^n)_{n \in \mathbb{N}}$ derived from solution $(\mathbf{y}^n, g^n)_{n \in \mathbb{N}}$ is not consistent with the continuous evolution. Even if the *Crank-Nicolson scheme* provides a second order in time variation of the total energy uniformly with respect to $\varepsilon > 0$ as indicated by (2.5) and (2.7), the variations of the discrete energy and discrete potential energy are not consistent.

Thus, several works have been devoted to modifications of the Crank-Nicolson scheme for (1.4) to obtain a result of uniform consistency with respect to $\varepsilon$. For instance, we mention the work of Brackbill, Forslund and Vu [5, 26] who first introduced an effective force into the equation on $\mathbf{v}$ in order to capture the $\nabla_{\mathbf{x}}^{\perp} b/b^2$
term in the limit $\epsilon \to 0$. We also refer to the recent work of Ricketson and Chacón [23], who proposed an alternative approach, which is expected to conserve energy within the Particle-In-Cell framework.

**2.1. The Brackbill-Forslund-Vu scheme.** The scheme developed by Brackbill, Forslund and Vu in [5,26] incorporates an effective force into the second equation of the Crank-Nicolson scheme (2.1). This additional force is designed to capture the correct drift $\nabla_x^\perp b/b^2$ when $\epsilon \to 0$. Obviously in the regime $\Delta t \ll \epsilon$, this force is expected to be significantly small, actually of order $O(\Delta t^2/\epsilon^4)$, see [5] for instance. More precisely, for a given time step $\Delta t > 0$, we define $t^n = n\Delta t$, for $n \in \mathbb{N}$ and $(\mathbf{x}_\epsilon^n, \mathbf{v}_\epsilon^n)$, an approximation of the solution $(\mathbf{x}_\epsilon, \mathbf{v}_\epsilon)$ to (1.4) at time $t^n$, through

$$ (2.9) \qquad \left\{ \begin{aligned} & \epsilon \frac{\mathbf{x}_\epsilon^{n+1} - \mathbf{x}_\epsilon^n}{\Delta t} = \mathbf{v}_\epsilon^{n+1/2}, \\ & \epsilon \frac{\mathbf{v}_\epsilon^{n+1} - \mathbf{v}_\epsilon^n}{\Delta t} = \mathbf{E}(\mathbf{x}_\epsilon^{n+1/2}) + \mathbf{F}_{\text{eff}}^{n+1/2} - b(\mathbf{x}_\epsilon^{n+1/2}) \frac{(\mathbf{v}_\epsilon^{n+1/2})^\perp}{\epsilon}, \\ & \mathbf{x}_\epsilon^0 = \mathbf{x}(0), \quad \mathbf{v}_\epsilon^0 = \mathbf{v}(0), \end{aligned} \right. $$

where the effective force $\mathbf{F}_{\text{eff}}$ is given by

$$ (2.10) \qquad \mathbf{F}_{\text{eff}}^{n+1/2} := -\eta^{n+1/2} \frac{\nabla_x b}{b}(\mathbf{x}_\epsilon^{n+1/2}), \quad \eta^{n+1/2} = \frac{1}{2} \left( \frac{\|\mathbf{v}_\epsilon^{n+1}\|^2}{2} + \frac{\|\mathbf{v}_\epsilon^n\|^2}{2} - \|\mathbf{v}_\epsilon^{n+1/2}\|^2 \right), $$

and we again use notation

$$ \mathbf{v}_\epsilon^{n+1/2} = \frac{\mathbf{v}_\epsilon^{n+1} + \mathbf{v}_\epsilon^n}{2}, \quad \mathbf{x}_\epsilon^{n+1/2} = \frac{\mathbf{x}_\epsilon^{n+1} + \mathbf{x}_\epsilon^n}{2}. $$

We now define the discrete kinetic energy $e_\epsilon^n = \frac{1}{2}\|\mathbf{v}_\epsilon^n\|^2$ and total energy as

$$ \mathcal{E}_\epsilon^n := e_\epsilon^n + \phi(\mathbf{x}_\epsilon^n), $$

in which $\phi$ is a given smooth function $\phi \in W^{3,\infty}(\mathbb{R}^2)$. Hence, we obtain the variation of the discrete kinetic energy by multiplying the second equation of (2.9) by $\mathbf{x}_\epsilon^{n+1/2}$, which gives

$$ (2.11) \qquad \frac{e_{\epsilon}^{n+1} - e_{\epsilon}^{n}}{\Delta t} = \mathbf{E}(\mathbf{x}_{\epsilon}^{n+1/2}) \cdot \frac{\mathbf{v}_{\epsilon}^{n+1/2}}{\epsilon} + \mathbf{F}_{\text{eff}}^{n+1/2} \cdot \frac{\mathbf{v}_{\epsilon}^{n+1/2}}{\epsilon}. $$

Then, applying a Taylor expansion to the potential $\phi$, it yields that

$$ (2.12) \qquad \begin{aligned} \frac{\mathcal{E}_{\epsilon}^{n+1} - \mathcal{E}_{\epsilon}^{n}}{\Delta t} &= \frac{e_{\epsilon}^{n+1} - e_{\epsilon}^{n}}{\Delta t} + \frac{\phi_{\epsilon}^{n+1} - \phi_{\epsilon}^{n}}{\Delta t}, \\ &= \mathbf{F}_{\text{eff}}^{n+1/2} \cdot \frac{\mathbf{x}_{\epsilon}^{n+1} - \mathbf{x}_{\epsilon}^{n}}{\Delta t} + \Delta t^2 \mathcal{O}\left(\left\|\frac{\mathbf{x}_{\epsilon}^{n+1} - \mathbf{x}_{\epsilon}^{n}}{\Delta t}\right\|^3\right). \end{aligned} $$

Using the definition of the effective force $\mathbf{F}_{\text{eff}}^{n+1/2}$, we have

$$ |\eta^{n+1/2}| = \frac{1}{8} \|\mathbf{v}^{n+1} - \mathbf{v}^n\|^2 = \left(\frac{\Delta t}{\epsilon}\right)^2 \mathcal{O}\left(\left\|\frac{\mathbf{x}_\epsilon^{n+1} - \mathbf{x}_\epsilon^n}{\Delta t}\right\|^3\right). $$

Then under our above assumptions and the further reasonable assumption that $\mathbf{x}_\epsilon^n$ is bounded, the evolution of the discrete energy obtained by (2.9) is much worse than the one (2.5) corresponding to the Crank-Nicolson scheme. Now, let us study the asymptotic behavior of (2.9) as $\epsilon$ tends to zero.

**Proposition 2.2** (Asymptotic behavior $\epsilon \to 0$ with a fixed $\Delta t$). Let $\phi \in W^{3,\infty}(\mathbb{R}^2)$, choose a sufficiently small fixed time step $\Delta t$ and a final time $T > 0$. We set $N_T = \lceil T/\Delta t \rceil$. Assume that the modified Crank-Nicolson scheme (2.9) defines a numerical approximation $(\mathbf{x}_\epsilon^n, \mathbf{v}_\epsilon^n)_{0 \le n \le N_T}$ satisfying

(i) for all $1 \le n \le N_T$, $\mathbf{x}_\epsilon^n$ is uniformly bounded with respect to $\epsilon > 0$;

(ii) in the limit $\epsilon \to 0$, $(\mathbf{x}_\epsilon^0, \frac{1}{2}\|\mathbf{v}_\epsilon^0\|^2)$ converges to some $(\mathbf{y}^0, g^0)$.

Then, we have
* for all $1 \le n \le N_T$, $(\mathbf{x}_\varepsilon^n, e_\varepsilon^n)$ converges to $(\mathbf{y}^n, g^n)$, as $\varepsilon \to 0$ with $e_\varepsilon^n = \frac{1}{2}\|\mathbf{v}_\varepsilon^n\|^2$ and the limit $(\mathbf{y}^n, g^n)_{1\le n\le N_T}$ solves

$$ (2.13) \qquad \begin{cases} \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}^{n+1/2}) + g^0 \frac{\nabla_\mathbf{y}^\perp b}{b^2}(\mathbf{y}^{n+1/2}), \\ \frac{g^{n+1} - g^n}{\Delta t} = 0. \end{cases} $$

* for all $1 \le n \le N_T$, the total energy $\mathcal{E}_\varepsilon^n = e_\varepsilon^n + \phi(\mathbf{x}_\varepsilon^n)$ converges to $\mathcal{E}_{gc}^n := g^n + \phi(\mathbf{y}^n)$ as $\varepsilon \to 0$, which satisfies

$$ (2.14) \qquad \frac{\mathcal{E}_{gc}^{n+1} - \mathcal{E}_{gc}^n}{\Delta t} = -g^0 \frac{\mathbf{E} \cdot \nabla_\mathbf{y}^\perp b}{b^2} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2); $$

* defining the discrete magnetic moment $\mu_\varepsilon^n = \frac{e_\varepsilon^n}{b(\mathbf{x}_\varepsilon^n)}$ converges to $\mu_{gc}^n := \frac{g^n}{b(\mathbf{y}^n)}$ as $\varepsilon \to 0$ such that

$$ (2.15) \qquad \frac{\mu_{gc}^{n+1} - \mu_{gc}^n}{\Delta t} = -g^0 \frac{\mathbf{E} \cdot \nabla_\mathbf{y}^\perp b}{b^3} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2). $$

where $\mathbf{y}^{n+1/2}$ is defined as $\mathbf{y}^{n+1/2} = (\mathbf{y}^{n+1} + \mathbf{y}^n)/2$.

*Proof.* The beginning of the proof of Proposition 2.1 applies word by word to the present case since it only uses the first equation of the scheme. In this way one arrives at a stage where one knows that $(\varepsilon^{-1}\mathbf{v}_\varepsilon^{n+1/2})_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon$ and a subsequence $(\mathbf{x}_\varepsilon^n, e_\varepsilon^n)$ converges to some $(\mathbf{y}^n, g^n)$ as $\varepsilon$ goes to zero, satisfying the second equation of (2.13).

Note that this also implies

$$ \lim_{\varepsilon \to 0} \eta^{n+1/2} = g^0 $$

so that when taking the limit $\varepsilon \to 0$ in the second equation of (2.9) one receives

$$ \lim_{\varepsilon \to 0} \frac{\mathbf{v}_\varepsilon^{n+1/2}}{\varepsilon} = -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}^{n+1/2}) + g^0 \frac{\nabla_\mathbf{y}^\perp b}{b^2}(\mathbf{y}^{n+1/2}). $$

This is sufficient to complete the derivation of (2.13). As in the proof of Proposition 2.1, we conclude the full convergence (and not only the convergence of a subsequence) from the fact that (2.13) defines a unique solution.

We then proceed as in the proof of Proposition 2.1 for the evolution of the total energy $\mathcal{E}_\varepsilon^n$ and the magnetic moment $\mu_\varepsilon^n$, and derive

$$ \frac{\mathcal{E}_{gc}^{n+1} - \mathcal{E}_{gc}^n}{\Delta t} = \nabla_{\mathbf{y}} \phi(\mathbf{y}^{n+1/2}) \cdot \left( \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} \right) + \mathcal{O}(\Delta t^2) = -g^0 \frac{\mathbf{E} \cdot \nabla_{\mathbf{y}}^\perp b}{b^2} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2) $$

and

$$ \frac{\mu_{gc}^{n+1} - \mu_{gc}^n}{\Delta t} = -g^0 \frac{(\mathbf{y}^{n+1} - \mathbf{y}^n)}{\Delta t} \cdot \frac{\nabla_{\mathbf{y}} b}{b^2} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2) = -g^0 \frac{\mathbf{E} \cdot \nabla_{\mathbf{y}}^\perp b}{b^3} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2). $$

□

Here it is worth mentioning that Proposition 2.2 shows that the scheme with the effective force (2.9) does not give a consistent approximation of slow variables $(\mathbf{x}_\varepsilon, e_\varepsilon)$ in the limit $\varepsilon \to 0$. Indeed, in the limit $\varepsilon \to 0$, the scheme (2.9) exactly preserves the kinetic energy over time while this quantity should vary according to the gradient of the magnetic field (grad **B** drift). Therefore, neither the discrete guiding center variable $\mathbf{y}^n$ nor the kinetic energy $g^n$ are consistent approximation of the guiding center system (1.8). As a consequence, the evolution of the discrete magnetic moment $\mu_{gc}$ is also not consistent with the continuous equation (1.10). Furthermore, the scheme (2.9) fails to preserve the second order accuracy with respect to $\Delta t$ of the total energy $\mathcal{E}_{gc}$, in contrast to the Crank-Nicolson scheme (2.1).

In order to overcome this drawback, an alternative numerical scheme, still based on the Crank-Nicolson method has been proposed recently by Ricketson and Chacón [23].
2.2. **The Ricketson-Chacón scheme.** The numerical scheme proposed by L. F. Ricketson and L. Chacón [23] still consists in adding a force term to capture the $\nabla_{\mathbf{x}}^{\perp}b/b^2$ drift in the asymptotic limit $\varepsilon \to 0$. More precisely, the force $\mathbf{F}_{\text{eff}}^{n+1/2}$ is chosen to be orthogonal to the velocity $\mathbf{v}_{\varepsilon}^{n+1/2}$, so that it does not interfere with the evolution of the discrete kinetic energy. Explicitly as in [23], set

$$ (2.16) \qquad \mathbf{F}_{\text{cons}}^{n+1/2} = \left( \text{Id} - \frac{\mathbf{v}_{\varepsilon}^{n+1/2} \otimes \mathbf{v}_{\varepsilon}^{n+1/2}}{\|\mathbf{v}_{\varepsilon}^{n+1/2}\|^2} \right) \mathbf{G}^{n+1/2}, $$

where $\mathbf{G}^{n+1/2}$ is given by

$$ (2.17) \qquad \mathbf{G}^{n+1/2} = \begin{cases} 2 \mathbf{F}_{\text{eff}}^{n+1/2} & \text{if } \|\mathbf{v}_{\varepsilon}^{n+1/2} - \mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2}\| \ge \|\mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2}\|, \\ \left( \frac{2}{\beta_{\varepsilon}^{n+1/2}} \hat{\mathbf{v}}_{\mathbf{E},\varepsilon}^{n+1/2} \otimes \hat{\mathbf{v}}_{\mathbf{E},\varepsilon}^{n+1/2} + \frac{\text{Id} - \hat{\mathbf{v}}_{\mathbf{E},\varepsilon}^{n+1/2} \otimes \hat{\mathbf{v}}_{\mathbf{E},\varepsilon}^{n+1/2}}{1 - \frac{\beta_{\varepsilon}^{n+1/2}}{2}} \right) \mathbf{F}_{\text{eff}}^{n+1/2} & \text{otherwise} \end{cases} $$

with

$$ \mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2} = -\frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}_{\varepsilon}^{n+1/2}), \quad \hat{\mathbf{v}}_{\mathbf{E},\varepsilon}^{n+1/2} = \frac{\mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2}}{\|\mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2}\|}, \quad \beta_{\varepsilon}^{n+1/2} = \frac{\|\mathbf{v}_{\varepsilon}^{n+1/2} - \mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2}\|^2}{\|\mathbf{v}_{\mathbf{E},\varepsilon}^{n+1/2}\|^2}, $$

whereas the effective force $\mathbf{F}_{\text{eff}}^{n+1/2}$ is given in (2.10) and $\text{Id}$ is the identity matrix. Then, the modified Crank-Nicolson scheme now becomes [23]

$$ (2.18) \qquad \left\{ \begin{aligned} \varepsilon \frac{\mathbf{x}_{\varepsilon}^{n+1} - \mathbf{x}_{\varepsilon}^{n}}{\Delta t} &= \mathbf{v}_{\varepsilon}^{n+1/2}, \\ \varepsilon \frac{\mathbf{v}_{\varepsilon}^{n+1} - \mathbf{v}_{\varepsilon}^{n}}{\Delta t} &= \mathbf{E}(\mathbf{x}_{\varepsilon}^{n+1/2}) + \mathbf{F}_{\text{cons}}^{n+1/2} - b(\mathbf{x}_{\varepsilon}^{n+1/2}) \frac{(\mathbf{v}_{\varepsilon}^{n+1/2})^{\perp}}{\varepsilon}, \\ \mathbf{x}_{\varepsilon}^{0} &= \mathbf{x}(0), \quad \mathbf{v}_{\varepsilon}^{0} = \mathbf{v}(0). \end{aligned} \right. $$

where again

$$ \mathbf{v}_{\varepsilon}^{n+1/2} = \frac{\mathbf{v}_{\varepsilon}^{n+1} + \mathbf{v}_{\varepsilon}^{n}}{2}, \quad \mathbf{x}_{\varepsilon}^{n+1/2} = \frac{\mathbf{x}_{\varepsilon}^{n+1} + \mathbf{x}_{\varepsilon}^{n}}{2}. $$

As we did previously, we define the kinetic energy $e_{\varepsilon}^{n} = \frac{1}{2}\|\mathbf{v}_{\varepsilon}^{n}\|^2$ and the discrete total energy as $\mathcal{E}_{\varepsilon}^{n} = e_{\varepsilon}^{n} + \phi(\mathbf{x}_{\varepsilon}^{n})$ for a given potential charge $\phi \in W^{3,\infty}(\mathbb{R}^2)$. Using that $\mathbf{F}_{\text{cons}}^{n+1/2}$ is orthogonal to $\mathbf{v}_{\varepsilon}^{n+1/2}$ thus to $\mathbf{x}_{\varepsilon}^{n+1} - \mathbf{x}_{\varepsilon}^{n}$, we recover the same evolution of the discrete total energy as the one for the Crank-Nicolson scheme: for all $\varepsilon > 0$,

$$ (2.19) \qquad \frac{\mathcal{E}_{\varepsilon}^{n+1} - \mathcal{E}_{\varepsilon}^{n}}{\Delta t} = \Delta t^2 \mathcal{O} \left( \left\| \frac{\mathbf{x}_{\varepsilon}^{n+1} - \mathbf{x}_{\varepsilon}^{n}}{\Delta t} \right\|^3 \right). $$

Note that strictly speaking, because of (2.16), when $\mathbf{v}_{\varepsilon}^{n+1/2}$ is zero an alternative for $\mathbf{F}_{\text{cons}}^{n+1/2}$ should be used. Now, let us investigate the asymptotic limit of the scheme (2.18) when $\varepsilon \to 0$ with a fixed $\Delta t$.

**Proposition 2.3** (Asymptotic behavior $\varepsilon \to 0$ with a fixed $\Delta t$). Let $\phi \in W^{3,\infty}(\mathbb{R}^2)$ such that $\nabla\phi$ is nowhere vanishing, choose a sufficiently small fixed time step $\Delta t$ and a final time $T > 0$. We set $N_T = \lfloor T/\Delta t \rfloor$. Assume that the modified Crank-Nicolson scheme (2.18) defines a numerical approximation $(\mathbf{x}_{\varepsilon}^{n}, \mathbf{v}_{\varepsilon}^{n})_{0 \le n \le N_T}$ satisfying

(i) for all $1 \le n \le N_T$, $\mathbf{x}_{\varepsilon}^{n}$ is uniformly bounded with respect to $\varepsilon > 0$;

(ii) in the limit $\varepsilon \to 0$, $(\mathbf{x}_{\varepsilon}^{0}, \frac{1}{2}\|\mathbf{v}_{\varepsilon}^{0}\|^2)$ converges to some $(\mathbf{y}^{0}, g^{0})$.

Then, we have

* for all $1 \le n \le N_T$, $(\mathbf{x}_{\varepsilon}^{n}, e_{\varepsilon}^{n})$ converges to $(\mathbf{y}^{n}, g^{n})$, as $\varepsilon \to 0$ with $e_{\varepsilon}^{n} = \frac{1}{2}\|\mathbf{v}_{\varepsilon}^{n}\|^2$ and the limit $(\mathbf{y}^{n}, g^{n})_{1 \le n \le N_T}$ solves

$$ (2.20) \qquad \left\{ \begin{aligned} \frac{\mathbf{y}^{n+1} - \mathbf{y}^{n}}{\Delta t} &= -\frac{\mathbf{E}^{\perp}}{b}(\mathbf{y}^{n+1/2}) + 2g^{0} \left( \frac{\mathbf{E} \cdot \nabla_{\mathbf{y}} b}{b^2 \|\mathbf{E}\|^2} \right) \mathbf{E}^{\perp}(\mathbf{y}^{n+1/2}), \\ \frac{g^{n+1} - g^{n}}{\Delta t} &= 0. \end{aligned} \right. $$
• for all $1 \le n \le N_T$, the total energy $\mathcal{E}_\varepsilon^n = e_\varepsilon^n + \phi(\mathbf{x}_\varepsilon^n)$ converges to $\mathcal{E}_{gc}^n := g^n + \phi(\mathbf{y}^n)$ as $\varepsilon \to 0$, which satisfies

$$ (2.21) \qquad \frac{\mathcal{E}_{gc}^{n+1} - \mathcal{E}_{gc}^n}{\Delta t} = \mathcal{O}(\Delta t^2), $$

• defining the discrete magnetic moment $\mu_\varepsilon^n = \frac{e_\varepsilon^n}{b(\mathbf{x}_\varepsilon^n)}$ converges to $\mu_{gc}^n := \frac{g^n}{b(\mathbf{y}^n)}$ as $\varepsilon \to 0$, which satisfies

$$ (2.22) \qquad \frac{\mu_{gc}^{n+1} - \mu_{gc}^n}{\Delta t} = -g^0 \frac{\mathbf{E}^{\perp} \cdot \nabla_y b}{b^2} (\mathbf{y}^{n+1/2}) - 2 (g^0)^2 \frac{(\mathbf{E}^{\perp} \cdot \nabla_y b) (\mathbf{E} \cdot \nabla_y b)}{b^4 \|\mathbf{E}\|^2} (\mathbf{y}^{n+1/2}) + \mathcal{O}(\Delta t^2). $$
where $\mathbf{y}^{n+1/2}$ is defined as $\mathbf{y}^{n+1/2} = (\mathbf{y}^{n+1} + \mathbf{y}^n)/2$.

**Proof.** The beginning of the proof of Proposition 2.2 applies verbatim to the present case, since it only utilizes the first equation of the scheme and the definition of $\mathbf{F}_{\text{eff}}^{n+1/2}$. In this way, one arrives at a stage where it is known that $(\varepsilon^{-1}\mathbf{v}_\varepsilon^{n+1/2})_{\varepsilon>0}$ is uniformly bounded with respect to $\varepsilon$ and a subsequence $(\mathbf{x}_\varepsilon^n, e_\varepsilon^n)$ converges to some $(\mathbf{y}^n, g^n)$ as $\varepsilon$ goes to zero, satisfying the second equation of (2.20) and

$$ \lim_{\varepsilon \to 0} \mathbf{F}_{\text{eff}}^{n+1/2} = -g^0 \frac{\nabla_y b}{b} (\mathbf{y}^{n+1/2}). $$

Note that, since by assumption $\mathbf{E}$ is nowhere vanishing, this also implies

$$ \lim_{\varepsilon \to 0} \beta^{n+1/2} = 1 $$

so that

$$ \lim_{\varepsilon \to 0} \mathbf{G}^{n+1/2} = 2 \lim_{\varepsilon \to 0} \mathbf{F}_{\text{eff}}^{n+1/2} = -2 g^0 \frac{\nabla_y b}{b} (\mathbf{y}^{n+1/2}). $$

This time the determination of the limit of $\varepsilon^{-1}\mathbf{v}_\varepsilon^{n+1/2}$ is much more complicated. By extracting further if necessary we may assume that it converges to some $\mathbf{u}^{n+1/2}$ and that $\mathbf{v}_\varepsilon^{n+1/2}/\|\mathbf{v}_\varepsilon^{n+1/2}\|$ converges to some $\hat{\mathbf{u}}^{n+1/2}$. By using that for any nonzero $\mathbf{z}$

$$ \mathrm{Id} - \frac{\mathbf{z} \otimes \mathbf{z}}{\|\mathbf{z}\|^2} = \frac{\mathbf{z}^{\perp} \otimes \mathbf{z}^{\perp}}{\|\mathbf{z}\|^2} $$

and taking a limit in the second equation of (2.18), we derive

$$ \hat{\mathbf{u}}^{n+1/2} \left( \|\mathbf{u}^{n+1/2}\| - 2g^0 \hat{\mathbf{u}}^{n+1/2} \cdot \frac{\nabla_y^\perp b}{b^2} (\mathbf{y}^{n+1/2}) \right) = - \frac{\mathbf{E}^\perp}{b} (\mathbf{y}^{n+1/2}) $$

Since $\mathbf{E}$ is non vanishing, this implies that $\hat{\mathbf{u}}^{n+1/2}$ is colinear to $\mathbf{E}^\perp$ and thus

$$ \mathbf{u}^{n+1/2} = -\frac{\mathbf{E}^{\perp}}{b}(\mathbf{y}^{n+1/2}) + 2g^0 \mathbf{E}^{\perp}(\mathbf{y}^{n+1/2}) \frac{\mathbf{E} \cdot \nabla_y b}{b^2 \|\mathbf{E}\|^2} (\mathbf{y}^{n+1/2}). $$

This completes the derivation of (2.20), which then may be used as before to upgrade the convergence to the full convergence.

The rest of the proof for the variations of the total discrete energy and the discrete adiabatic invariant is then analogous to the ones in the proof of Propositions 2.1 and 2.2. $\square$

This latter Proposition shows that the modified scheme (2.18) does not provide a consistent asymptotic limit when $\varepsilon$ goes to zero and $\Delta t$ is fixed. Actually, in [23], the authors propose an adaptive time step procedure to overcome this drawback.

**Remark 2.4.** It is worth mentioning that the Crank-Nicolson is well suited to design an approximation preserving the total energy. In particular, we refer to [24] where the following modified electric field

$$ (2.23) \qquad \tilde{\mathbf{E}}(\mathbf{x}_\varepsilon^{n+1/2}) = \frac{\phi(\mathbf{x}_\varepsilon^{n+1}) - \phi(\mathbf{x}_\varepsilon^n)}{(\mathbf{x}_\varepsilon^{n+1} - \mathbf{x}_\varepsilon^n) \cdot \mathbf{E}(\mathbf{x}_\varepsilon^{n+1/2})} \mathbf{E}(\mathbf{x}_\varepsilon^{n+1/2}). $$

is applied ensuring exact preservation of the total energy, that is $\mathcal{E}_\varepsilon^n = \mathcal{E}_\varepsilon^0$, for all $n \in \mathbb{N}$. However, this exact preservation does not help to provide a consistent approximation as $\varepsilon \to 0$ since the discrete kinetic energy is not uniformly consistent with respect to $\varepsilon$.
### 3. A MODIFIED CRANK-NICOLSON SCHEME WITH AN ADDITIONAL VARIABLE

We now propose a new numerical scheme based on the Crank-Nicolson method (2.1), designed to be asymptotically consistent with the guiding center model (1.8) as $\epsilon \to 0$. To this aim, we apply the strategy developed in [12] which consists in solving an augmented system incorporating the discrete kinetic energy $(e_\epsilon^n)_{n \in \mathbb{N}}$ into the discrete system. Furthermore, as in the previous work, we add an effective force to capture the drift $\nabla_{\mathbf{x}}^{\perp} b/b^2$ in the limit $\epsilon \to 0$. More precisely, we reformulate the system (1.4) for $(\mathbf{x}_\epsilon, \mathbf{v}_\epsilon)$ in an equivalent manner for the new unknowns $(\mathbf{x}_\epsilon, \mathbf{w}_\epsilon, e_\epsilon)$ as

$$ (3.1) \qquad \left\{ \begin{aligned} \epsilon \frac{d\mathbf{x}_\epsilon}{dt} &= \mathbf{w}_\epsilon, \\ \epsilon \frac{d\mathbf{e}_\epsilon}{dt} &= \mathbf{E}(\mathbf{x}_\epsilon) \cdot \mathbf{w}_\epsilon, \\ \epsilon \frac{d\mathbf{w}_\epsilon}{dt} &= \mathbf{E}(\mathbf{x}_\epsilon) - \chi(\mathbf{w}_\epsilon, e_\epsilon) \frac{\nabla_{\mathbf{x}} b}{b}(\mathbf{x}_\epsilon) - b(\mathbf{x}_\epsilon) \frac{\mathbf{w}_\epsilon^\perp}{\epsilon}, \\ \mathbf{x}_\epsilon^0 &= \mathbf{x}(0), \quad \mathbf{w}_\epsilon^0 = \mathbf{v}(0), \quad e_\epsilon^0 = \frac{1}{2} \|\mathbf{v}(0)\|^2, \end{aligned} \right. $$

where $\chi$ is chosen as

$$ \chi(\mathbf{w}, e) = \max \left( e - \frac{1}{2} \|\mathbf{w}\|^2, 0 \right), \quad \forall (\mathbf{w}, e) \in \mathbb{R}^2 \times \mathbb{R}^+. $$

Then, we discretize this system applying a classical Crank-Nicolson scheme to $(\mathbf{x}_\epsilon^n, \mathbf{w}_\epsilon^n, e_\epsilon^n)$,

$$ (3.2) \qquad \left\{ \begin{aligned} \epsilon \frac{\mathbf{x}_\epsilon^{n+1} - \mathbf{x}_\epsilon^n}{\Delta t} &= \mathbf{w}_\epsilon^{n+1/2}, \\ \epsilon \frac{e_\epsilon^{n+1} - e_\epsilon^n}{\Delta t} &= \mathbf{E}(\mathbf{x}_\epsilon^{n+1/2}) \cdot \mathbf{w}_\epsilon^{n+1/2}, \\ \epsilon \frac{\mathbf{w}_\epsilon^{n+1} - \mathbf{w}_\epsilon^n}{\Delta t} &= \mathbf{E}(\mathbf{x}_\epsilon^{n+1/2}) - \chi(\mathbf{w}_\epsilon^{n+1/2}, e_\epsilon^{n+1/2}) \frac{\nabla_{\mathbf{x}} b}{b}(\mathbf{x}_\epsilon^{n+1/2}) - b(\mathbf{x}_\epsilon^{n+1/2}) \frac{(\mathbf{w}_\epsilon^{n+1/2})^\perp}{\epsilon}, \\ \mathbf{x}_\epsilon^0 &= \mathbf{x}_\epsilon(0), \quad \mathbf{w}_\epsilon^0 = \mathbf{v}_\epsilon(0), \quad e_\epsilon^0 = \frac{1}{2} \|\mathbf{v}_\epsilon(0)\|^2, \end{aligned} \right. $$

in which

$$ \mathbf{w}_\epsilon^{n+1/2} = \frac{\mathbf{w}_\epsilon^{n+1} + \mathbf{w}_\epsilon^n}{2}, \quad \mathbf{x}_\epsilon^{n+1/2} = \frac{\mathbf{x}_\epsilon^{n+1} + \mathbf{x}_\epsilon^n}{2}. $$

Computing an approximation at time $t^{n+1}$ requires the numerical solution of a nonlinear system. Here, we employ a straightforward fixed-point iteration scheme, where the variable $\mathbf{x}_\epsilon^{n+1/2}$, depending on $\mathbf{x}_\epsilon^{n+1}$, is frozen at each step. This approach reduces the problem to solving a two-dimensional linear system for the velocity variable $\mathbf{w}_\epsilon^{n+1}$ of each particle. The convergence tolerance is set to $10^{-10}$. Although alternative approaches based on Newton's method could be utilized, the current fixed-point scheme converges rapidly, typically requiring fewer than five iterations.

At each time step, the velocity $(\mathbf{v}_\epsilon^n)_{n \in \mathbb{N}}$ is given by

$$ \mathbf{v}_\epsilon^n = \sqrt{2e_\epsilon^n} \frac{\mathbf{w}_\epsilon^n}{\|\mathbf{w}_\epsilon^n\|}. $$

As for the original Crank-Nicolson scheme (2.5), the variation of the discrete total energy obeys

$$ (3.3) \qquad \frac{\mathcal{E}_\epsilon^{n+1} - \mathcal{E}_\epsilon^n}{\Delta t} = \frac{e_\epsilon^{n+1} - e_\epsilon^n}{\Delta t} + \frac{\phi_\epsilon^{n+1} - \phi_\epsilon^n}{\Delta t} = \Delta t^2 \mathcal{O}\left(\left\|\frac{\mathbf{x}_\epsilon^{n+1} - \mathbf{x}_\epsilon^n}{\Delta t}\right\|^3\right). $$

Therefore, under the reasonable assumption that $\mathbf{x}_\epsilon^n$ is bounded, the variation of the total discrete energy is of order $\Delta t^2$. Now, let us investigate the asymptotic behavior of the scheme (3.2) as $\epsilon$ goes to zero.

**Proposition 3.1** (Consistency in the limit $\epsilon \to 0$ for a fixed $\Delta t$). Let $\phi \in W^{3,\infty}(\mathbb{R}^2)$, choose an a priori bound $M$ and a final time $T > 0$ then a sufficiently small fixed time step $\Delta t$. We set $N_T = \lfloor T/\Delta t \rfloor$. Assume that the Crank-Nicolson scheme (3.2) defines a numerical approximation $(\mathbf{x}_\epsilon^n, \mathbf{w}_\epsilon^n, e_\epsilon^n)_{0 \le n \le N_T}$ satisfying

(i) for all $1 \le n \le N_T$, $\mathbf{x}_\epsilon^n$ is uniformly bounded with respect to $\epsilon > 0$;

(ii) in the limit $\epsilon \to 0$, $(\mathbf{x}_\epsilon^0, e_\epsilon^0)$ converges to some $(\mathbf{y}^0, g^0)$ such that $g^0 \le M$.
Then we have

• for all $1 \le n \le N_T$, $(\mathbf{x}_\epsilon^n, e_\epsilon^n)$ converges to $(\mathbf{y}^n, g^n)$, as $\epsilon \to 0$ and the limit $(\mathbf{y}^n, g^n)_{1\le n\le N_T}$ solves

$$ (3.4) \qquad \begin{cases} \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}^{n+1/2}) + g^{n+1/2} \frac{\nabla_\mathbf{y}^\perp b}{b^2}(\mathbf{y}^{n+1/2}), \\ \frac{g^{n+1} - g^n}{\Delta t} = g^{n+1/2} \frac{\mathbf{E} \cdot \nabla_\mathbf{y}^\perp b}{b^2}(\mathbf{y}^{n+1/2}). \end{cases} $$

• for all $1 \le n \le N_T$, the total energy $\mathcal{E}_\epsilon^n = e_\epsilon^n + \phi(\mathbf{x}_\epsilon^n)$ converges to $\mathcal{E}_{gc}^n := g^n + \phi(\mathbf{y}^n)$ as $\epsilon \to 0$, which satisfies

$$ (3.5) \qquad \frac{\mathcal{E}_{gc}^{n+1} - \mathcal{E}_{gc}^n}{\Delta t} = \mathcal{O}(\Delta t^2), $$

• defining the discrete magnetic moment $\mu_\epsilon^n = \frac{e_\epsilon^n}{b(\mathbf{x}_\epsilon^n)}$ converges to $\mu_{gc}^n := \frac{g^n}{b(\mathbf{y}^n)}$ as $\epsilon \to 0$ such that

$$ (3.6) \qquad \frac{\mu_{gc}^{n+1} - \mu_{gc}^n}{\Delta t} = \mathcal{O}(\Delta t^2) $$

where $(\mathbf{y}^{n+1/2}, g^{n+1/2})$ is defined as $\mathbf{y}^{n+1/2} = (\mathbf{y}^{n+1} + \mathbf{y}^n)/2$ and $g^{n+1/2} = (g^{n+1} + g^n)/2$.

**Proof.** We again follow the lines of the proof of Proposition 2.1. To begin with, from the first line of (3.2) we deduce that $(\epsilon^{-1}\mathbf{w}_\epsilon^{n+1/2})_{\epsilon>0}$ is uniformly bounded with respect to $\epsilon$. Combined with the second line of (3.2) this implies that each $(e_\epsilon^{n+1/2})_{\epsilon>0}$ is also uniformly bounded with respect to $\epsilon$.

Therefore, up to a subsequence $(\mathbf{x}_\epsilon^n, e_\epsilon^n)$ converges to $(\mathbf{y}^n, g^n)$ as $\epsilon$ goes to zero. One readily deduces that

$$ \lim_{\epsilon \to 0} \chi(\mathbf{w}_\epsilon^{n+1/2}, e_\epsilon^{n+1/2}) = g^{n+1/2}, $$

so that from the third equation of (3.2) stems

$$ \lim_{\epsilon \to 0} \frac{\mathbf{w}_\epsilon^{n+1/2}}{\epsilon} = -\frac{\mathbf{E}^\perp}{b}(\mathbf{y}^{n+1/2}) + g^{n+1/2} \frac{\nabla_\mathbf{y}^\perp b}{b^2}(\mathbf{y}^{n+1/2}). $$

Inserting the latter in the first and second equations of (3.2) completes the derivation of (3.4). Then, again, the convergence is upgraded from the convergence of a subsequence to full convergence by the uniqueness of solutions to the limiting system, (3.4). Note that the system (3.4) is more nonlinear from previously derived asymptotic systems, which is why the required constraint on $\Delta t$ depends here on a priori bound on $g^0$ and $T$.

The rest of the proof, on the variations of the total discrete energy and the discrete magnetic moment, is omitted as completely analogous to the corresponding one of Proposition 2.1. $\square$

**Proposition 3.1** shows that the modified scheme (3.2) is consistent uniformly with respect to $\epsilon$ and allows to recover a consistent approximation of the guiding center system (1.8). Moreover, the new scheme also preserves the second order accuracy for the total energy and the magnetic moment.

Let us also stress that a straightforward adaptation of the conservation trick (2.23) from [24] provides a genuinely energy conserving version of the present scheme. Yet our numerical simulations, not reported here, show no further significant improvement so that we have decided not to delve further into this direction.

It is worth noting that this scheme does not guarantee the nonnegativity of the variable $e_\epsilon^n$ for all $n \in \mathbb{N}$. If such a situation arises, one possible remedy is to update $e_\epsilon^n$ as

$$ e_\epsilon^n = \frac{1}{2} \|\mathbf{w}_\epsilon^n\|^2. $$

However, in the upcoming numerical simulations, this scenario never occurs.

## 4. NUMERICAL SIMULATIONS

In this section, we provide examples of numerical computations to validate and compare the different time discretization schemes introduced in the previous sections. We first consider the motion of a single particle under the effect of a given electromagnetic field. It allows us to illustrate the theoretical results in the limit $\epsilon \to 0$ of the numerical schemes and their accuracy for multi-scale problems.

Then we consider the Vlasov-Poisson system with an external non uniform magnetic field. We apply a classical Particle-In-Cell method with the time discretization technique based on the Crank-Nicolson
scheme (3.2) to describe the diocotron instability in a disk and also the stability of vortices in a D-shape domain.

4.1. **One single particle motion.** We first investigate the motion of an individual particle in a given electromagnetic field. We consider the electric field $\mathbf{E} = -\nabla_x \phi$ where the potential $\phi$ is given by

$$ (4.1) \qquad \phi(\mathbf{x}) = \frac{x_2^2}{2}, $$

while the external magnetic field is

$$ (4.2) \qquad b(\mathbf{x}) = 1 + \|\mathbf{x}\|^2. $$

The initial condition is chosen as $\mathbf{x}^0 = (2,2)$, $\mathbf{v}^0 = (3,3)$ and the final time $T = 1$. On the one hand, we compute reference solutions $(\mathbf{x}_\varepsilon, \mathbf{w}_\varepsilon, e_\varepsilon)_{\varepsilon>0}$ and $(\mathbf{y}, g)$ to the stiff initial value problem (1.4) and to the asymptotic problem (1.8) thanks to an explicit fourth-order Runge-Kutta scheme using a small time step chosen according to the size of order $O(\varepsilon^2)$ for the initial system. On the other hand, for various time steps $\Delta t$, independent of $\varepsilon$, we compute approximate solutions $(\mathbf{x}_{\varepsilon,\Delta t}, \mathbf{w}_{\varepsilon,\Delta t}, e_{\varepsilon,\Delta t})$ using the modified Crank-Nicolson scheme (3.2) and also compare the results with those obtained using (2.9) proposed in [5, 26], and (2.18) described in [7, 23]. For completeness, we also compare our results with those obtained using an IMEX2L for the augmented system (1.6) developed in [11, 12]. To evaluate the accuracy, the numerical error is measured as

$$ \left\{ \begin{aligned} \|\mathbf{x}_{\varepsilon,\Delta t} - \mathbf{x}_{\varepsilon}\| &:= \frac{\Delta t}{T} \sum_{n=0}^{N_T} \|\mathbf{x}_{\varepsilon,\Delta t}^n - \mathbf{x}_{\varepsilon}(t^n)\|, \\ \|\mathbf{x}_{\varepsilon,\Delta t} - \mathbf{y}\| &:= \frac{\Delta t}{T} \sum_{n=0}^{N_T} \|\mathbf{x}_{\varepsilon,\Delta t}^n - \mathbf{y}(t^n)\|, \\ \|e_{\varepsilon,\Delta t} - g\| &:= \frac{\Delta t}{T} \sum_{n=0}^{N_T} |e_{\varepsilon,\Delta t}^n - g(t^n)|. \end{aligned} \right. $$

In Figure 4.1, we present the numerical error on $\|\mathbf{x}_{\varepsilon,\Delta t} - \mathbf{x}_{\varepsilon}\|$ expressed with respect to $\varepsilon$ in log-log scale for various time steps $\Delta t \in \{10^{-5}, \dots, 10^{-1}\}$. When $\varepsilon \ge 10^{-1}$, we observe the expected second order accuracy of the different schemes. However, as $\varepsilon$ becomes smaller, the time steps are too large and the numerical error for both schemes (2.9) and (2.18) increases. In contrast, the numerical error associated to the modified Crank-Nicolson scheme (3.2) and the IMEX2L [12] decreases with respect to $\varepsilon$. This behavior is typical of an asymptotic preserving scheme for which the error becomes of order $\varepsilon$ when the time step $\Delta t$ is sufficiently large. It is worth mentioning that the behavior of the errors for (3.2) and the IMEX2L differs significantly as $\varepsilon \ll 1$. Indeed, even with a large time step the modified Crank-Nicolson scheme (3.2) remains so accurate that the error of order $\Delta t^2$ is negligible compared to the error with respect to $\varepsilon$. This phenomenon can also be observed in Figures 4.2 and 4.3, where we report the errors compared to the reference solution of the asymptotic model $\|\mathbf{x}_{\varepsilon,\Delta t} - \mathbf{y}\|$ and $\|e_{\varepsilon,\Delta t} - g\|$. Clearly the schemes (2.9) and (2.18) do not capture a consistent approximation $(\mathbf{y}, g)$ to the asymptotic solutions (1.8) as $\varepsilon \to 0$. These numerical experiments illustrate the lack of consistency proven in Propositions 2.2 and 2.3. In contrast, both schemes (3.2) and the IMEX2L, based on the approximation of the augmented system (1.6), successfully capture the limit with the correct convergence rate with respect to $\varepsilon$ (slope of order one). The advantage of the scheme (3.2) is that, when $\varepsilon \ll 1$, the amplitude of the numerical error is much smaller than the one corresponding to other schemes.

To illustrate this point, we also present the space trajectories corresponding to $\varepsilon = 0.01$ and $\Delta t = 0.1$ for large time simulations (with $T = 30$). In Figure 4.4, we observe that the particle trajectory forms a circular motion under the effect of both drifts $\mathbf{E}^{\perp}/b$ and $\nabla_x^{\perp}b/b^2$, which appear explicitly in the asymptotic model (1.8). The scheme (3.2) provides an approximation close to the reference solution, while other schemes fail to get the correct position of the particle, since they do not capture correctly the drift $\nabla_x^{\perp}b/b^2$ in the limit $\varepsilon \to 0$.

Furthermore, from the results presented in Figure 4.5, we observe that the IMEX2L and (3.2) schemes accurately track the variations in kinetic and potential energy over a long period, unlike (2.9) and (2.18).

Finally, the time evolution of the magnetic moment approximation is presented and compared with a reference solution in Figure 4.6. Let us emphasize that only the modified scheme (3.2) accurately describes the amplitude of the fast oscillations of $\mu_\varepsilon$ compared to the other schemes. Indeed, even if the time step
<figure><img src="image_2.png" /></figure>

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

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

<figure><img src="image_5.png" /><figcaption>FIGURE 4.1. **One single particle motion:** Numerical errors of discrete solution $\mathbf{x}_{\varepsilon, \Delta t}$, approximated by several schemes: (2.9), (2.18), IMEX2L and (3.2), with reference solution $\mathbf{x}_{\varepsilon}$ of (1.4) for various $\varepsilon > 0$ and $\Delta t > 0$.</figcaption></figure>

$\Delta t$ is much larger than the fastest time scale of order $\varepsilon^2$, the quantity $\mu_{\varepsilon, \Delta t}$ oscillates with the correct amplitude of order $\varepsilon^2$.

These numerical results on the single particle motion, clearly show that the modified Crank-Nicolson scheme (3.2) is much more accurate than the schemes (2.9), (2.18) or IMEX2L, even when $\varepsilon \ll 1$, with a fixed time step $\Delta t$ independent of $\varepsilon$. These elementary numerical simulations confirm the ability of the modified Crank-Nicolson scheme (3.2) to capture the evolution of the "slow" variables $(\mathbf{x}_{\varepsilon}, e_{\varepsilon})$ uniformly with respect to $\varepsilon$ by essentially transitioning automatically to the guiding center motion (1.8) when it is not able anymore to follow the fast oscillations of the initial system.

Let us remind that the scheme proposed in [23] uses both an effective force and an adaptive time-stepping procedure as $\varepsilon$ becomes small. Our goal here is to deepen the understanding of the underlying issues, particularly by distinguishing the effects of the adaptive time step procedure from those of using an effective force. In particular, we illustrate the theoretical results obtained in Sections 2 and 3 and demonstrate that the inclusion of a slow variable allows us to dispense with the adaptive time step procedure, enabling the use of an arbitrarily large time step independent of $\varepsilon$. Applying the adaptive time step procedure allows
<figure><img src="image_6.png" /></figure>

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

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

<figure><img src="image_9.png" /><figcaption>FIGURE 4.2. **One single particle motion:** Numerical errors of discrete solution $x_{\epsilon, \Delta t}$, approximated by several schemes: (2.9), (2.18), IMEX2L and (3.2), with guiding center solution y of (1.8) for various $\epsilon > 0$ and $\Delta t > 0$.</figcaption></figure>

significantly to decrease the error when $\epsilon$ goes to zero, but it requires a small time step depending on $\epsilon$, which affects the computational cost.

To end this section dedicated to comparing the various numerical schemes, Table 1 presents the computational cost of the IMEX2L scheme, which does not require an iterative method, alongside that of the modified Crank-Nicolson scheme (3.2). Since the number of iterations does not exceed 5 at each time step for all $\Delta t$ and $\epsilon$, it is observed that the computational time for scheme (3.2) is more than twice that of the IMEX2L scheme. However, the error associated with (3.2) is significantly smaller, especially when $\epsilon \ll 1$.

Let us conclude this section with an important remark about strongly oscillating fields.

**Remark 4.1.** The scheme (3.2) and the strategy proposed in [13] can also be applied to the case where the electric field is highly oscillatory, that is, $||\partial_t E||/||E|| = O(1/\epsilon)$. In such a situation, the asymptotic limit remains given by (1.8), and the scheme (3.2) can be applied directly. Indeed, we performed numerical
<figure><img src="image_10.png" /><figcaption>FIGURE 4.3. One single particle motion: Numerical errors of discrete solution $e_{\epsilon, \Delta t}$, approximated by several schemes: (2.9), (2.18), IMEX2L and (3.2), with guiding center solution $g$ of (1.8) for various $\epsilon > 0$ and $\Delta t > 0$.</figcaption></figure>

<table><caption>TABLE 1. Comparison of the computational time (microsecond) of the scheme (3.2) with IMEX2L, with final time $T = 1$ for various $\epsilon > 0$ and $\Delta t > 0$.</caption><thead><tr><th></th><th>Δt = 10<sup>-1</sup></th><th>Δt = 10<sup>-2</sup></th><th>Δt = 10<sup>-3</sup></th><th>Δt = 10<sup>-4</sup></th></tr></thead><tbody><tr><td>ε = 10<sup>-1</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>Scheme (3.2)</td><td>5</td><td>34</td><td>314</td><td>2572</td></tr><tr><td>IMEX2L</td><td>2</td><td>17</td><td>176</td><td>1663</td></tr><tr><td>ε = 10<sup>-2</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>Scheme (3.2)</td><td>5</td><td>65</td><td>289</td><td>3399</td></tr><tr><td>IMEX2L</td><td>2</td><td>17</td><td>177</td><td>1803</td></tr><tr><td>ε = 10<sup>-3</sup></td><td></td><td></td><td></td><td></td></tr><tr><td>Scheme (3.2)</td><td>5</td><td>31</td><td>265</td><td>2617</td></tr><tr><td>IMEX2L</td><td>2</td><td>17</td><td>183</td><td>1656</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>FIGURE 4.4. One single particle motion: Trajectory of particle approximated by several schemes: (2.9), (2.18), IMEX2L and (3.2) with $\epsilon = 0.01$, $\Delta t = 0.1$ and final time $T = 30s$.</figcaption></figure>

simulations (not presented in this paper) with a potential

$$ \phi(t, \mathbf{x}) = \frac{1}{2} \cos \left( \frac{t}{\epsilon} \right) \|\mathbf{x}\|^2 $$

and obtained the same error curves as those presented in this section.

## 4.2. Vlasov-Poisson system.
We now consider the Vlasov-Poisson system (1.3) on a domain $\Omega \subset \mathbb{R}^2$, where $\Omega$ is given either by a disk or a D shape domain (see Figure 4.7).

Assuming that the density is concentrated far from the boundary, we choose to remove particles which are located outside the physical domain, this may induce a lack of mass conservation. However, thanks to the strong confinement of the magnetic field, we do not observe this situation in the present simulations. For the potential $\phi_{\epsilon}$, the Poisson equation is solved with homogeneous Dirichlet boundary conditions using a classical five points second order finite difference method with ghost points to take into account the effect of the boundary conditions [28].
<figure><img src="image_12.png" /><figcaption>FIGURE 4.5. One single particle motion: Variation of energy approximated by several schemes: (2.9), (2.18), IMEX2L and (3.2) with $\epsilon = 0.01$, $\Delta t = 0.1$ and final time $T = 30s$.</figcaption></figure>

For each time $t \ge 0$, we may define the total energy $\mathcal{E}_\epsilon(t)$ as

$$ (4.3) \qquad \mathcal{E}_{\epsilon}(t) := \mathbf{K}_{\epsilon}(t) + \mathbf{U}_{\epsilon}(t), $$

where the kinetic energy $\mathbf{K}_\epsilon(t)$ and the potential energy $\mathbf{U}_\epsilon(t)$ are given by

$$ \left\{ \begin{array}{l} \mathbf{K}_{\epsilon}(t) := \frac{1}{2} \int_{\Omega} \int_{\mathbb{R}^2} f_{\epsilon}(t, \mathbf{x}, \mathbf{v}) \|\mathbf{v}\|^2 d\mathbf{v} d\mathbf{x}, \\ \mathbf{U}_{\epsilon}(t) := \frac{1}{2} \int_{\Omega} \|\mathbf{E}_{\epsilon}(t, \mathbf{x})\|^2 d\mathbf{x}. \end{array} \right. $$

Assuming that the distribution function is compactly supported in the open set $\Omega$, the total energy $\mathcal{E}_\epsilon(t)$ is conserved for all time $t \ge 0$.
<figure><img src="image_13.png" /><figcaption>FIGURE 4.6. One single particle motion: Evolution of magnetic moment approximated by several schemes: (2.9), (2.18), IMEX2L and (3.2) with $\epsilon = 0.01$, $\Delta t = 0.1$ and final time $T = 30s$.</figcaption></figure>

We also define the magnetic moment for the Vlasov-Poisson system (1.3) given by

$$ (4.4) \qquad \mu_{\epsilon}(t) := \frac{1}{2} \iint_{\Omega} \int_{\mathbb{R}^2} f_{\epsilon}(t, \mathbf{x}, \mathbf{v}) \frac{\|\mathbf{v}\|^2}{b(\mathbf{x})} d\mathbf{v} d\mathbf{x}, \quad t \ge 0 $$

and expect that $\mu_{\epsilon}(t)$ is an invariant in time in the asymptotic limit $\epsilon \to 0$ for the limit model (1.8).

For this section, we performed numerical experiments using the modified Crank-Nicolson scheme (3.2) to approximate the particles trajectory corresponding to the Vlasov equation. Despite the fact that the modified Crank-Nicolson scheme (3.2) does not conserves exactly the total energy, we expect that its variations are of order $O(\Delta t^2)$ even when $\epsilon$ tends to zero. Furthermore, the modified Crank-Nicolson scheme (3.2) should capture correctly the asymptotic limit $\epsilon \to 0$, as it has been shown for the single particle motion.

### 4.2.1. Diocotron instability.
We first consider Vlasov-Poisson system (1.3) set in a disk $\Omega = D(0, 12)$ centered at the origin with a radius $R = 12$ as in Figure 4.7. Here, the Particle-In-Cell method is implemented
<figure><img src="image_14.png" /><figcaption>FIGURE 4.7. Disk domain (left) and D-Shape domain (right) embedded in uniform Cartersian grid.</figcaption></figure>

with approximately 100 particles per cell on a uniform grid of the square $(-12, 12)^2$ with $\Delta x = 0.1$. The simulation starts with a Maxwellian distribution in velocity, whose macroscopic density is a perturbed uniform distribution in an annulus. More precisely, we choose

$$f(0, \mathbf{x}, \mathbf{v}) = \frac{\rho_0(\mathbf{x})}{2\pi} \exp\left(-\frac{\|\mathbf{v}\|^2}{2}\right),$$

where $\rho_0$ is given by

$$\rho_0(\mathbf{x}) = \begin{cases} n_0(1 + \alpha \cos(7\theta)), & \text{for } 6 \le \|\mathbf{x}\| \le 7, \\ 0, & \text{else,} \end{cases}$$

in which $n_0 = 0.25$, $\alpha = 0.001$, and the angle $\theta$ is defined as $\theta = \arctan(x_2/x_1)$ with $\mathbf{x} = (x_1, x_2)$. Moreover, we consider strong external magnetic field which is given by

$$b(t, \mathbf{x}) = \frac{20}{\sqrt{400 - \|\mathbf{x}\|^2}}.$$

Since $b$ is not homogeneous, even in the asymptotic regime the kinetic and potential parts of the total energy are not preserved separately, but the total energy corresponding to the Vlasov-Poisson system is still preserved. Figure 4.8 shows that all these features are captured satisfactorily by the modified Crank-Nicolson scheme (3.2) even on long time evolutions with a large time step $\Delta t = 0.1$ and small $\varepsilon = 10^{-2}$. On the one hand, the variations of the total energy have an amplitude of order $10^{-4}$, which is satisfying compared to the physical variations of the potential and kinetic energy of order $2 \times 10^{-2}$. On the other hand, the quantity $\mu_\varepsilon$ also varies around $10^{-4}$, which corresponds to the scale of $\varepsilon^2 = 10^{-4}$. This phenomenon has already been observed for the single particle motion and will be discussed below.

In Figure 4.9, we visualize the corresponding dynamics by presenting several snapshots of the macroscopic density at some specific times $t = 0, 50, 100$ and $150$. The numerical results obtained with our PIC methods are in good agreements with those obtained with a finite difference scheme [17] for the asymptotic model (1.11).

Furthermore, we performed simulations for various $\varepsilon \in \{10^{-2}, 5 \cdot 10^{-2}, 10^{-1}\}$, shown in Figure 4.10. On the left hand side, we report the total energy variations, which are theoretically of order $\Delta t^2$. However, when $\varepsilon$ is small these variations decrease so that the total energy is well preserved in the limit $\varepsilon \to 0$. On the right hand side, we present the variations of the adiabatic invariant $\mu_\varepsilon$, which is not preserved by the solution to the Vlasov-Poisson system but only by the asymptotic model (1.11). Here, we notice that this quantity oscillates with an amplitude of order $\varepsilon^2$. Surprisingly, even with a large time step $\Delta t$, the numerical scheme (3.2) is able to capture the correct amplitude. This can also be observed on the variations
<figure><img src="image_15.png" /><figcaption>FIGURE 4.8. Diocotron instability: Time evolution of the variations of the relative potential $U_{\epsilon}$ and kinetic energy $K_{\epsilon}$ (left) and magnetic moment $\mu_{\epsilon}$ (right) with $\epsilon = 10^{-2}$ and $(\Delta t, \Delta x) = (0.1, 0.1)$, using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

of both the potential and the kinetic energy in Figures 4.8. The slow variations definitively correspond to the effect of the drifts $E^{\perp}/b$ and $\nabla_{\mathbf{x}}^{\perp}b/b^2$, but the fast oscillations and their amplitudes are more intricate. In order to verify that these small oscillations are not a numerical artefact, we compute a reference solution with a small time step for $\epsilon = 10^{-1}$ and $10^{-2}$ and compare these results with those obtained from (3.2) with $\Delta t = 0.1$. Figure 4.11 clearly indicates that with such a time step, the scheme (3.2) described accurately the amplitude of these fast oscillations. However, since $\Delta t$ is much larger than the oscillation period, the modified scheme can not describe the physical frequency.

4.2.2. Vortex interaction. Finally we consider Vlasov-Poisson system (1.3) in a D-Shaped domain $\Omega \subset \mathbb{R}^2$, described by a mapping from polar coordinates $(r, \theta)$ to Cartesian coordinates $\mathbf{x} = (x_1, x_2)$ given by

$$ \begin{cases} x_1 = a + r \cos(\theta + \arcsin(0.416) \sin(\theta)), \\ x_2 = b + 1.66 r \sin(\theta), \end{cases} $$

centered at the origin $(a, b) = (0, 0)$ where $0 < r \le 10$ and $0 \le \theta \le 2\pi$ as in Figure 4.7. Here, we consider the Particle-In-Cell method with approximatively 100 particles per cell on a uniform grid of the rectangle $(-11, 11) \times (-17, 17)$ with a space discretization $\Delta \mathbf{x} = 0.1$. We choose the initial distribution function such that

$$ (4.5) \qquad f(0, \mathbf{x}, \mathbf{v}) = \frac{5}{8\pi^2} \left[ \exp\left(-\frac{\|\mathbf{x} - \mathbf{x}_0\|^2}{2}\right) + \exp\left(-\frac{\|\mathbf{x} + \mathbf{x}_0\|^2}{2}\right) \right] \exp\left(-\frac{\|\mathbf{v}\|^2}{2}\right), $$

with $\mathbf{x}_0 = (1.5, -1.5)$. Moreover, we consider a non homogeneous external magnetic field such as

$$ b(t, \mathbf{x}) = \frac{20}{\sqrt{400 - x_1^2 - x_2^2}}. $$

As expected for such a configuration, since $b$ is not homogeneous, even in the asymptotic regime the kinetic and potential parts of the total energy are not preserved separately, but the total energy corresponding to the Vlasov-Poisson system is still preserved. In addition, the quantity $\mu_{\epsilon}$ is an invariant for the guiding center model but oscillates with a high frequency for the Vlasov-Poisson system with a small amplitude. Indeed, Figure 4.12 shows that all these features are again captured by the scheme (3.2) even on long time evolutions with a large time step. Furthermore, in Figure 4.13, we visualize the corresponding dynamics by presenting several snapshots of the time evolution of the macroscopic charge density for $\epsilon = 10^{-2}$ at time $t = 0, 80, 160, 240, 320$ and $400$. Since $\epsilon \ll 1$, the conservation of $e/b(\mathbf{x})$ offers coercivity jointly in $(\mathbf{x}, e)$ allowing to confine the density in the D-shape domain. Such a confinement is indeed observed, jointly with the expected eventual merging of two initial vortices in a relatively short time. We also
<figure><img src="image_16.png" /><figcaption>FIGURE 4.9. **Diocotron instability:** Macroscopic density evolution at some specific time $T = 0, 50, 100, 150$ with $\epsilon = 10^{-2}$ and $(\Delta t, \Delta x) = (0.1, 0.1)$, using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

observe small filaments at low density, which generate a "halo" propagating into the domain as already observed in [?, 17]. Finally, from Figure 4.14, similarly to diocotron instability experiment in Section 4.2.1, the relative variations of the total energy $\mathcal{E}_\epsilon$ and the adiabatic invariant $\mu_\epsilon$ show the ability of preserving these parameters $(\mathcal{E}_\epsilon, \mu_\epsilon)$ in the limit $\epsilon \to 0$.

In summary, the modified Crank-Nicolson scheme (3.2) integrated into Particle-In-Cell method shows consistency and stability of long time simulations with a coarse time step, even for small $\epsilon \ll 1$. Moreover, the solutions $(\mathcal{E}_\epsilon, \mu_\epsilon)$ of these numerical tests preserve the structure of the limit system $\epsilon \to 0$ and provide accurately the amplitude of these variations.

## 5. CONCLUSION

In this paper, we propose a modified Crank-Nicolson time discretization technique for Particle-In-Cell simulations. Our approach guarantees the accuracy and stability of small-scale variables, even when the magnetic field amplitude becomes large, thus correctly capturing their long-term behavior, including in cases of inhomogeneous magnetic fields and coarse time grids. Comparison with previous contributions on Crank-Nicolson and semi-implicit schemes demonstrates the effectiveness of our approach, with accuracy improving by several orders of magnitude. Even for large-time simulations, the resulting numerical schemes provide acceptable accuracy on physical invariants (total energy for all $\epsilon$, magnetic moment when $\epsilon \ll 1$), while fast scales are automatically filtered when the time step is large compared to $\epsilon^2$.
<figure><img src="image_17.png" /><figcaption>FIGURE 4.10. **Diocotron instability:** Time evolution of the variations of the relative total energy $\mathcal{E}_{\varepsilon}$ (left) and magnetic moment $\mu_{\varepsilon}$ (right) with different $\varepsilon = 10^{-1}, 5 \cdot 10^{-2}$ and $10^{-2}$ with $(\Delta t, \Delta x) = (0.1, 0.1)$, using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 4.11. **Diocotron instability:** Time evolution of the variation of the potential energy $U_{\varepsilon}$ in a short time interval with $\varepsilon = 0.1$ (left) and $\varepsilon = 0.01$ (right) using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

As a theoretical validation, we have proven that, under certain stability assumptions on the numerical approximations, the slow part of the approximation converges, when $\varepsilon \to 0$, to the solution of a limit scheme of asymptotic evolution, preserving the initial order of precision. Let us notice that here we chose to use a Crank-Nicolson scheme, which is very popular in computational physics; however, our strategy can be easily applied to higher-order schemes [11] (such as IMEX or multi-step methods).

The next step involves extending this strategy to 3D particle simulations, incorporating curvature effects and addressing both the parallel and orthogonal directions relative to the magnetic field. These aspects give rise to more complex phenomena; however, by considering the slow variables and the conservation of energy, we are able to apply the same approach effectively [14]. A complete analytical study (including proofs of stability), in the spirit of [16], is currently under investigation [15].
<figure><img src="image_19.png" /><figcaption>FIGURE 4.12. Vortex interaction: Time evolution of the variations of the relative potential $U_{\epsilon}$ and kinetic energy $K_{\epsilon}$ (left) and magnetic moment $\mu_{\epsilon}$ (right) with $\epsilon = 10^{-2}$ with $(\Delta t, \Delta x) = (0.1, 0.1)$, using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

## APPENDIX A. FORMAL ASYMPTOTIC BEHAVIOR FOR A GIVEN ELECTROMAGNETIC FIELD

For the convenience of the reader, we provide here the main algebraic manipulations of equations (1.4) leading to the guiding-center system (1.8). The system is well-known and no details of the derivation is needed in the rest of the paper but consistency with (1.8) is crucially used in our evaluations of numerical schemes.

To begin with, in order to eliminate $\epsilon^{-1}\mathbf{v}_{\epsilon}$ from (1.6), we observe that equation (1.7) may be replaced with

$$ (A.1) \qquad \epsilon \frac{d}{dt} \left( \frac{\mathbf{v}_{\epsilon}^{\perp}}{b(\mathbf{x}_{\epsilon})} \right) = \frac{\mathbf{E}^{\perp}(\mathbf{x}_{\epsilon})}{b(\mathbf{x}_{\epsilon})} - \left( \frac{\nabla_{\mathbf{x}} b(\mathbf{x}_{\epsilon})}{b^2(\mathbf{x}_{\epsilon})} \cdot \mathbf{v}_{\epsilon} \right) \mathbf{v}_{\epsilon}^{\perp} + \frac{\mathbf{v}_{\epsilon}}{\epsilon}. $$

In order to characterize the asymptotic dynamics of the slow variables $(\mathbf{x}_{\epsilon}, e_{\epsilon})$ when $\epsilon \to 0$, we notice that their equations depend linearly on $\mathbf{v}_{\epsilon}$, hence as in [13] we write for all $t \in \mathbb{R}^{+}$ and any linear operator $\mathbf{L}(t)$

$$ \epsilon \frac{d}{dt} \left( \mathbf{L}(t) \frac{\mathbf{v}_{\epsilon}^{\perp}}{b(\mathbf{x}_{\epsilon})} \right) = \epsilon \frac{d\mathbf{L}(t)}{dt} \frac{\mathbf{v}_{\epsilon}^{\perp}}{b(\mathbf{x}_{\epsilon})} + \mathbf{L}(t) \left( \frac{\mathbf{E}^{\perp}(\mathbf{x}_{\epsilon})}{b(\mathbf{x}_{\epsilon})} - \left( \frac{\nabla_{\mathbf{x}} b(\mathbf{x}_{\epsilon})}{b^2(\mathbf{x}_{\epsilon})} \cdot \mathbf{v}_{\epsilon} \right) \mathbf{v}_{\epsilon}^{\perp} + \frac{\mathbf{v}_{\epsilon}}{\epsilon} \right). $$

Therefore, applying the latter to $\mathbf{L}_{\mathbf{x}}(t) : \mathbf{u} \in \mathbb{R}^{2} \mapsto \mathbf{u} \in \mathbb{R}^{2}$ and $\mathbf{L}_{\epsilon}(t) : \mathbf{u} \in \mathbb{R}^{2} \mapsto \mathbf{E}(\mathbf{x}_{\epsilon}) \cdot \mathbf{u} \in \mathbb{R}$ and inserting the outcome into the system of slow variables (1.6), we obtain

$$ (A.2) \qquad \left\{ \begin{array}{l} \displaystyle \frac{d}{dt} \left( \mathbf{x}_{\epsilon} - \epsilon \frac{\mathbf{v}_{\epsilon}^{\perp}}{b(\mathbf{x}_{\epsilon})} \right) = - \frac{\mathbf{E}^{\perp}(\mathbf{x}_{\epsilon})}{b(\mathbf{x}_{\epsilon})} + \left( \frac{\nabla_{\mathbf{x}} b(\mathbf{x}_{\epsilon})}{b^2(\mathbf{x}_{\epsilon})} \cdot \mathbf{v}_{\epsilon} \right) \mathbf{v}_{\epsilon}^{\perp}, \\ \\ \displaystyle \frac{d}{dt} \left( e_{\epsilon} - \epsilon \mathbf{E}(\mathbf{x}_{\epsilon}) \cdot \frac{\mathbf{v}_{\epsilon}^{\perp}}{b(\mathbf{x}_{\epsilon})} \right) = - \left( \mathbf{v}_{\epsilon} \nabla_{\mathbf{x}} \right) \left( \frac{\mathbf{E}}{b} \right) (\mathbf{x}_{\epsilon}) \cdot \mathbf{v}_{\epsilon}^{\perp}. \end{array} \right. $$

This latter system may replace (1.6) and is not stiff with respect to $\epsilon \ll 1$, but it also suggests the introduction of new variables, as the guiding center variable $\mathbf{x}_{\epsilon} - \epsilon \mathbf{v}_{\epsilon}^{\perp}/b(\mathbf{x}_{\epsilon})$. However, these new equations still involve at leading order quadratic terms in $\mathbf{v}_{\epsilon}$ on the right hand side. Hence, following [13], we turn
<figure><img src="image_20.png" /><figcaption>FIGURE 4.13. Vortex interaction: Macroscopic density evolution at some specific time T with ε = 10⁻² with (Δt, Δx) = (0.1, 0.1), using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

our attention to bilinear operators $\mathbf{A}(t)$ and derive

$$ \begin{aligned} \text{(A.3)} \quad \varepsilon^2 \frac{d}{dt} \left( \mathbf{A}(t) \left( \mathbf{v}_{\varepsilon}, \frac{\mathbf{v}_{\varepsilon}^{\perp}}{b(\mathbf{x}_{\varepsilon})} \right) \right) &= \varepsilon^2 \frac{d\mathbf{A}(t)}{dt} \left( \mathbf{v}_{\varepsilon}, \frac{\mathbf{v}_{\varepsilon}^{\perp}}{b(\mathbf{x}_{\varepsilon})} \right) + \varepsilon \mathbf{A}(t) \left( \frac{\mathbf{E}(\mathbf{x}_{\varepsilon})}{b(\mathbf{x}_{\varepsilon})}, \mathbf{v}_{\varepsilon}^{\perp} \right) \\ &\quad + \varepsilon \mathbf{A}(t) \left( \mathbf{v}_{\varepsilon}, \frac{\mathbf{E}^{\perp}(\mathbf{x}_{\varepsilon})}{b(\mathbf{x}_{\varepsilon})} \right) - \varepsilon \left( \frac{\nabla_{\mathbf{x}} b(\mathbf{x}_{\varepsilon})}{b^2(\mathbf{x}_{\varepsilon})} \cdot \mathbf{v}_{\varepsilon} \right) \mathbf{A}(t) (\mathbf{v}_{\varepsilon}, \mathbf{v}_{\varepsilon}^{\perp}) \\ &\quad - \mathbf{A}(t) (\mathbf{v}_{\varepsilon}^{\perp}, \mathbf{v}_{\varepsilon}^{26}) + \mathbf{A}(t) (\mathbf{v}_{\varepsilon}, \mathbf{v}_{\varepsilon}). \end{aligned} $$
<figure><img src="image_21.png" /><figcaption>FIGURE 4.14. **Vortex interaction:** Time evolution of the variations of the relative total energy $\mathcal{E}_\varepsilon$ (left) and magnetic moment $\mu_\varepsilon$ (right) with different $\varepsilon = 10^{-1}, 5 \cdot 10^{-2}$ and $10^{-2}$ and $(\Delta t, \Delta x) = (0.1, 0.1)$, using the modified Crank-Nicolson scheme (3.2).</figcaption></figure>

Since $(\mathbf{v}_\varepsilon, \mathbf{v}_\varepsilon^\perp)$ is an orthogonal basis of $\mathbb{R}^2$ (when $\mathbf{v}_\varepsilon$ is non zero), introducing the operator Tr one observes that

$$ \|\mathbf{v}_\varepsilon\|^2 \mathrm{Tr}(\mathbf{A}(t)) = \mathbf{A}(t)(\mathbf{v}_\varepsilon^\perp, \mathbf{v}_\varepsilon^\perp) + \mathbf{A}(t)(\mathbf{v}_\varepsilon, \mathbf{v}_\varepsilon). $$

Therefore, one may reformulate (A.3) as

$$ \mathbf{A}(t)(\mathbf{v}_\varepsilon, \mathbf{v}_\varepsilon) = \frac{1}{2} \|\mathbf{v}_\varepsilon\|^2 \mathrm{Tr}(\mathbf{A}(t)) + \varepsilon^2 \frac{d\kappa_{\mathbf{A}}}{dt}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon) + \varepsilon \eta_{\mathbf{A}}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon), $$

where $\kappa_{\mathbf{A}}$ and $\eta_{\mathbf{A}}$ are given by

$$ \left\{ \begin{aligned} \kappa_{\mathbf{A}}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon) &= \frac{1}{2} \mathbf{A}(t) \left( \mathbf{v}_\varepsilon, \frac{\mathbf{v}_\varepsilon^\perp}{b(\mathbf{x}_\varepsilon)} \right), \\ \eta_{\mathbf{A}}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon) &= \frac{\nabla_{\mathbf{x}} b(\mathbf{x}_\varepsilon)}{2b^2(\mathbf{x}_\varepsilon)} \cdot \mathbf{v}_\varepsilon \mathbf{A}(t) \left( \mathbf{v}_\varepsilon, \mathbf{v}_\varepsilon^\perp \right) - \frac{\varepsilon}{2} \frac{d\mathbf{A}(t)}{dt} \left( \mathbf{v}_\varepsilon, \frac{\mathbf{v}_\varepsilon^\perp}{b(\mathbf{x}_\varepsilon)} \right) \\ &\quad - \frac{1}{2} \left( \mathbf{A}(t) \left( \frac{\mathbf{E}(\mathbf{x}_\varepsilon)}{b(\mathbf{x}_\varepsilon)}, \mathbf{v}_\varepsilon^\perp \right) + \mathbf{A}(t) \left( \mathbf{v}_\varepsilon, \frac{\mathbf{E}^\perp(\mathbf{x}_\varepsilon)}{b(\mathbf{x}_\varepsilon)} \right) \right). \end{aligned} \right. $$

To apply the latter to the bilinear maps $\mathbf{A}_x(t) : (\mathbf{u}_1, \mathbf{u}_2) \in \mathbb{R}^2 \times \mathbb{R}^2 \mapsto \mathbf{u}_2^\perp \frac{\nabla_{\mathbf{x}} b(\mathbf{x}_\varepsilon)}{b^2(\mathbf{x}_\varepsilon)} \cdot \mathbf{u}_1 \in \mathbb{R}^2$ and $\mathbf{A}_e(t) : (\mathbf{u}_1, \mathbf{u}_2) \in \mathbb{R}^2 \times \mathbb{R}^2 \mapsto (\mathbf{u}_1 \cdot \nabla_{\mathbf{x}}) \left( \frac{\mathbf{E}(\mathbf{x}_\varepsilon)}{b(\mathbf{x}_\varepsilon)} \right) \cdot \mathbf{u}_2^\perp \in \mathbb{R}$, we compute

$$ \mathrm{Tr}(\mathbf{A}_x(t)) = \frac{\nabla_{\mathbf{x}}^\perp b(\mathbf{x}_\varepsilon)}{b^2(\mathbf{x}_\varepsilon)} \quad \text{and} \quad \mathrm{Tr}(\mathbf{A}_e(t)) = \mathrm{div}_{\mathbf{x}} \left( -\frac{\mathbf{E}^\perp}{b}(\mathbf{x}_\varepsilon) \right). $$

Then we may replace (A.2) with the new

$$ \text{(A.4)} \quad \left\{ \begin{aligned} \frac{d}{dt} \left( \mathbf{x}_\varepsilon - \varepsilon \frac{\mathbf{v}_\varepsilon^\perp}{b(\mathbf{x}_\varepsilon)} - \varepsilon^2 \kappa_{\mathbf{x}_\varepsilon}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon) \right) &= -\frac{\mathbf{E}^\perp(\mathbf{x}_\varepsilon)}{b(\mathbf{x}_\varepsilon)} + e_\varepsilon \frac{\nabla_{\mathbf{x}}^\perp b}{b^2(\mathbf{x}_\varepsilon)}(\mathbf{x}_\varepsilon) + \varepsilon \eta_{\mathbf{x}_\varepsilon}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon), \\ \frac{d}{dt} \left( e_\varepsilon - \varepsilon \mathbf{E}(\mathbf{x}_\varepsilon) \cdot \frac{\mathbf{v}_\varepsilon^\perp}{b(\mathbf{x}_\varepsilon)} + \varepsilon^2 \kappa_{e_\varepsilon}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon) \right) &= e_\varepsilon \mathrm{div}_{\mathbf{x}} \left( \frac{\mathbf{E}^\perp}{b}(\mathbf{x}_\varepsilon) \right) - \varepsilon \eta_{e_\varepsilon}(t, \mathbf{x}_\varepsilon, \mathbf{v}_\varepsilon), \end{aligned} \right. $$

coupled with (1.7) for $\mathbf{v}_\varepsilon$ and $\kappa_\alpha$ and $\eta_\alpha$, for $\alpha \in \{\mathbf{x}_\varepsilon, e_\varepsilon\}$ are short-hand for $\kappa_{\mathbf{A}_\alpha}$ and $\kappa_{\mathbf{A}_\alpha}$.

This last formulation easily allows to characterize the asymptotic limit as $\varepsilon \to 0$ of the slow variables $(\mathbf{x}_\varepsilon, e_\varepsilon)$ provided one already knows that the fast variable $\mathbf{v}_\varepsilon$ remains bounded. Note in particular that $\mathbf{E}$ is curl-free thus $\mathbf{E}^\perp$ is divergence-free.
REFERENCES

[1] P. M. Bellan. *Fundamentals of plasma physics*. Cambridge University Press, 2008.

[2] C. K. Birdsall and A. B. Langdon. Plasma physics via computer simulation. *Series in plasma physics*. Taylor and Francis, New York, 2005. Originally published: New York ; London : McGraw-Hill, 1985.

[3] J. P. Boris. Relativistic Plasma Simulation - Optimization of a Hybrid code. *Proceeding of Fourth Conference on Numerical Simulations of Plasmas*, 1970.

[4] S. Boscarino, F. Filbet, and G. Russo. High order semi-implicit schemes for time dependent partial differential equations. *Journal of Scientific Computing*, 68:975-1001, 2016.

[5] J. U. Brackbill and D. W. Forslund. Simulation of Low-Frequency, Electromagnetic Phenomena in Plasmas. *Multiple Time Scales*, Academic Press, pages 271-310, 1985.

[6] M. Campos Pinto and F. Charles. From particle methods to forward-backward lagrangian schemes. *The SMAI Journal of computational mathematics*, 4:121-150, 2016.

[7] G. Chen and L. Chacón. An implicit, conservative and asymptotic-preserving electrostatic particle-in-cell algorithm for arbitrarily magnetized plasmas in uniform magnetic fields. *Journal of Computational Physics*, 487:112-160, 2023.

[8] N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field. *Journal of Computational Physics*, 346:172-190, 2017.

[9] P. Degond and F. Deluzet. Asymptotic-Preserving methods and multiscale models for plasma physics. *Journal of Computational Physics*, 336:429-457, 2017.

[10] P. Degond and F. Filbet. On the asymptotic limit of the three dimensional Vlasov-Poisson system for large magnetic field: Formal derivation. *Journal of Statistical Physics*, 165:765-784, 2016.

[11] F. Filbet and L. Rodrigues. Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field. *SIAM Journal on Numerical Analysis*, 54(2), 2016.

[12] F. Filbet and L. Rodrigues. Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas. *SIAM Journal on Numerical Analysis*, 55(5), 2017.

[13] F. Filbet and L. Rodrigues. Asymptotics of the three-dimensional Vlasov equation in the large magnetic field limit. *Journal of Ecole Polytechnique - Mathematics*, 7:1009-1067, 2020.

[14] F. Filbet and L. Rodrigues. Asymptotically preserving particle methods for strongly magnetized plasmas in a torus. *Journal of Computational Physics*, 480, 2023.

[15] F. Filbet, L. Rodrigues, and K. Trinh. Convergence analysis of a crank-nicolson scheme for strongly magnetized plasma. 2025.

[16] F. Filbet, L. M. Rodrigues, and H. Zakerzadeh. Convergence analysis of asymptotic preserving schemes for strongly magnetized plasmas. *Numerische Mathematik*, 149(3):549-593, 2021.

[17] F. Filbet and T. Xiong. Conservative discontinuous galerkin/hermite spectral method for the vlasov-poisson system. *Communications on Applied Mathematics and Computation*, 4(1):34-59, 2022.

[18] F. Filbet and C. Yang. Numerical simulations to the Vlasov-Poisson system with a strong magnetic field. arXiv:1805.10888, 2018.

[19] T. Genoni, R. Clark, and D. Welch. A fast implicit algorithm for highly magnetized charged particle motion. *The Open Plasma Physics Journal*, 3:36-41, 2010.

[20] R. Hazeltine and F. Waelbrock. The framework of plasma physics. CRC Press, 2018.

[21] R. D. Hazeltine and J. D. Meiss. Plasma confinement. *Courier Corporation*, 2003.

[22] K. Miyamoto. *Plasma physics and controlled nuclear fusion*, volume 38 of Springer Series on Atomic, Optical, and Plasma Physics. Springer-Verlag Berlin-Heidelberg, 2006.

[23] S. Parker and C. Birdsall. Numerical error in electron orbits with large $\omega_{ce}\delta t$. *Journal of Computational Physics*, 97:91-102, 1991.

[24] L. Ricketson and L. Chacón. An energy-conserving and asymptotic-preserving charged-particle orbit implicit time integrator for arbitrary electromagnetic fields. *Journal of Computational Physics*, 418, 2020.

[25] J. Simo, N. Tarnow, and K. Wong. Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics. *Computer Methods in Applied Mechanics and Engineering*, 100:63-116, 1992.

[26] A.-T. Vu. *Mathematical study of the Vlasov-Poisson equations with intense magnetic field*. PhD thesis, Aix-Marseille Université, 2023.

[27] H. X. Vu and J. U. Brackbill. Accurate numerical solution of charged particle motion in a magnetic field. *Journal of computational physics*, 116:384-387, 1995.

[28] B. Wang and X. Zhao. Error estimates of some splitting schemes for charged-particle dynamics under strong magnetic field. *SIAM Journal on Numerical Analysis*, 59:2075-2105, 2021.

[29] C. Yang and F. Filbet. Conservative and non-conservative methods based on hermite weighted essentially non-oscillatory reconstruction for vlasov equations. *Journal of Computational Physics*, 279:18-36, 2014.

UNIVERSITÉ DE, TOULOUSE, INSTITUT DE MATHÉMATIQUES DE TOULOUSE, FRANCE
*Email address: francis.filbet@math.univ-toulouse.fr*

UNIVERSITÉ DE RENNES, CNRS, IRMAR - UMR 6625, F-35000 RENNES, FRANCE
*Email address: luis-miguel.rodrigues@univ-rennes.fr*

UNIVERSITÉ DE RENNES, CNRS, IRMAR - UMR 6625, F-35000 RENNES, FRANCE
*Email address: kim-han.trinh@univ-rennes.fr*