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

HAL
open science

# Convergence analysis of a Crank-Nicolson scheme for strongly magnetized plasmas

Francis Filbet, L Miguel Rodrigues, Kim Han Trinh

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

## To cite this version:

Francis Filbet, L Miguel Rodrigues, Kim Han Trinh. Convergence analysis of a Crank-Nicolson scheme for strongly magnetized plasmas. 2025. hal-05394091

HAL Id: hal-05394091

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

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
CONVERGENCE ANALYSIS OF A CRANK-NICOLSON SCHEME
FOR STRONGLY MAGNETIZED PLASMAS

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

ABSTRACT. The present paper is devoted to the convergence analysis of an asymptotic preserving particle scheme designed to serve as a particle pusher in a Particle-In-Cell (PIC) method for the Vlasov equation with a strong inhomogeneous magnetic field. The asymptotic preserving scheme that we study removes classical strong restrictive stability constraints on discretization steps while capturing the large-scale dynamics, even when the discretization is too coarse to capture fastest scales. Our error bounds are explicit regarding the discretization and stiffness parameters and match sharply numerical tests. The present analysis is expected to be representative of the general analysis of a class of schemes, developed by the authors, conceived as implicit-explicit schemes on augmented formulations.

**Keywords:** Asymptotic preserving schemes; High-order time discretization; Strong magnetic field; Inhomogeneous magnetic field; Particle methods; Vlasov system.

**AMS Subject Classifications:** 65M75, 76X05, 65L04, 65M15, 35Q83, 82D10.

CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. Continuous analysis</td><td>7</td></tr><tr><td>3. Discrete uniform bounds</td><td>12</td></tr><tr><td>4. Asymptotic analysis</td><td>15</td></tr><tr><td>5. Direct convergence analysis</td><td>18</td></tr><tr><td>6. Proof of the main theorem, Theorem 1.1</td><td>25</td></tr><tr><td>Appendix A. Implicit solving</td><td>26</td></tr><tr><td>References</td><td>27</td></tr></table>

1. INTRODUCTION

The numerical analysis of the present paper is motivated by the numerical computations of the dynamics of plasmas where the charged particles evolve under an electrostatic and intense confining magnetic field. This configuration is typical of a tokamak plasma [1, 17] where the magnetic field is used to confine particles inside the core of the device. Kinetic models, based on a mesoscopic description of the various particles constituting a plasma, and coupled to Maxwell's equations for the computation of the electromagnetic fields, are very precise approaches for the study of such thermonuclear fusion plasmas. Incidentally we point out that magnetized plasmas are also encountered in a wide variety of astrophysical situations, and similar kinetic/Maxwell models, thus also our associated analysis, are relevant for those too.

Date: December 2, 2025.

Research of L.M.R. was partially supported by the ANR Project HEAD ANR-24-CE40-3260 and the Institut Universitaire de France.
Our present analysis contributes to the convergence analysis of numerical approximations obtained by particle methods (see [2]), which consist in approximating the particle distribution function by a finite number of macro-particles. The trajectories of these particles are determined from the characteristic curves associated to the underlying Vlasov equation; see (1.1) below for the characteristic differential equations studied here. From the presence of a strong magnetic field stems a stiff /dynamics where high-frequency oscillations co-exist with slower dynamics. Insisting on capturing the full stiff dynamics introduces a stringent restriction on discretization meshes. Over the last decade a wealth of energy has been poured into designing schemes relaxing this constraint. Let us simply refer the reader to [9] for a detailed discussion and comparison of various particle pushers for PIC methods and to the introduction of [8] for an account of other families of methods. In contrast, much fewer effort has been devoted to provide a complete mathematical analysis of the designed schemes. We mention [3, 12, 11, 19, 18, 20] as representative notable exceptions. It is worth pointing out that, in contrast with our present analysis, inspired by [10], all of those analyses fit in the frame of and use tools from the classical theory of geometric numerical integration, as described in [13].

Before delving into concrete description and analysis let us replace the scheme we propose and completely analyze here into the family of schemes [5, 6, 8, 9] it belongs to. The starting point is the hope that if one aims at simply capturing accurately the slow part of the dynamics the stiffness step constraint could be relaxed and the observation that in simplest cases this could be achieved through suitable semi-implicit (IMEX) particles pushers. In [5] this was carried out to capture the spatial dynamics in uniform magnetic fields. A complete analysis of this situation is provided in [10]. A genuinely new idea was needed to extend the approach to realistic situations when one cannot focus simply on the spatial dynamics. The solution proposed in [6], and explained further below, is to apply the approach of [5] to an extended system. In [6] this is achieved for inhomogeneous magnetic fields with a fixed direction. In [8] the strategy is adapted to more realistic geometries and in [9] it is shown how to improve conservation properties. Our present goal is to provide a complete analysis for schemes of the kind designed in [9], thus the present contribution is thought to be to [9] what [10] is to [5]. The scheme we completely analyze is very close to the ones introduced and tested in [9] but we propose a new scheme so as to minimize the technicalities of the analysis. We restrict here to proving error bounds for the computation of characteristics but we refer to [10] for a discussion on how this may be converted into errors on particle distribution functions.

Considering time-independent external electromagnetic fields with a fixed-direction magnetic field and an electric field deriving from an electric potential and focusing on the transverse dynamics on a time scale adapted to the strength of the magnetic field leads to the following set of characteristic differential system

$$ (1.1) \qquad \left\{ \begin{array}{l} \epsilon \frac{d\mathbf{x}}{dt} = \mathbf{v}, \\ \epsilon \frac{d\mathbf{v}}{dt} = \mathbf{E}(\mathbf{x}) - b(\mathbf{x}) \frac{\mathbf{v}^{\perp}}{\epsilon}, \end{array} \right. $$

where the transverse space-velocity unknown ($\mathbf{x}, \mathbf{v}$) is valued in $\mathbf{R}^2 \times \mathbf{R}^2$, $\epsilon > 0$ encodes the typical strength of the magnetic field, $b$ tracks its spatial variation and is lower bounded away from zero

$$ b(\mathbf{x}) \ge b_0 > 0, $$

$\mathbf{E}$ denotes the electric field and is assumed to take the form

$$ \mathbf{E}(\mathbf{x}) = -\nabla_{\mathbf{x}}\phi(\mathbf{x}) $$

for some electric potential $\phi$. Both $b$ and $\phi$ and all their derivatives are assumed to be smooth and bounded. The symbol $\perp$ denotes the direct orthogonal rotation, $\mathbf{w}^{\perp} = (-w_2, w_1)$. System (1.1)
is meant to be solved on a time interval $[0, T]$ and completed by an initial condition $\mathbf{x}(0) = \mathbf{x}^0$, $\mathbf{v}(0) = \mathbf{v}^0$. Note that solutions to (1.1) also satisfy the energy conservation

$$ \frac{d}{dt} \left( \frac{1}{2} \|\mathbf{v}\|^2 + \phi(\mathbf{x}) \right) = 0. $$

Since the arguments are to be generalized at the discrete level we provide in Section 2 a proof of the quite well-known fact, stated here as Proposition 2.3, that for any solution $(\mathbf{x}, \mathbf{v})$ to (1.1), $(\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2)$ remains $\varepsilon$-close to a solution $(\mathbf{y}, g)$ to

$$ (1.2) \qquad \left\{ \begin{aligned} \frac{d\mathbf{y}}{dt} &= -\frac{\mathbf{E}^{\perp}}{b}(\mathbf{y}) - g \nabla_{\mathbf{x}}^{\perp} \left(\frac{1}{b}\right)(\mathbf{y}), \\ \frac{dg}{dt} &= -\frac{d\phi(\mathbf{y})}{dt}. \end{aligned} \right. $$

When doing so we follow closely the approach of [7]. The variables $(\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2)$ are precisely the slow variables that we want to keep computing accurately even in situations when $\mathbf{v}$ is oscillating too fast to be reasonably computed. Before moving to the discrete level, let us stress that something that is difficult to enforce at the discrete level is that at the continuous level, as $\varepsilon$ goes to zero, the velocity $\mathbf{v}$ converges weakly to zero whereas the kinetic energy $\frac{1}{2}\|\mathbf{v}\|^2$ converges to a nonzero $g$ compatible with energy conservation, both being needed in the derivation of (1.2).

The accurate computation of $(\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2)$ shall be enforced through the requirement that the discretization should be compatible with $\varepsilon$ asymptotics. This is the asymptotic-preserving paradigm. For recent overviews on the latter we refer to [4, 16]. To achieve this, the relatively simple strategy initiated in [6] is to compute solutions to (1.2) as special solutions of a larger system where $(\mathbf{v}, \frac{1}{2}\|\mathbf{v}\|^2)$ are artificially separated, say as $(\mathbf{w}, e)$ so that numerical schemes are allowed to damp $\mathbf{w}$ to zero while letting $e$ converge to a suitable $g$ as $\varepsilon \to 0$. Note that the point is that the point is that the continuous time evolution preserves for the augmented vector $(x, e, \mathbf{w})$ preserves the relation $e = \frac{1}{2}\|\mathbf{w}\|^2$ whereas the corresponding discretized dynamics does not. Thus picking a numerical scheme for (1.1) along the strategy of [6, 8, 9] consists in choosing

* an augmented formulation for $(\mathbf{x}, e, \mathbf{w})$, reducing to (1.1) with $\mathbf{v} = \mathbf{w}$ when $e = \frac{1}{2}\|\mathbf{w}\|^2$

* a discretization of the augmented formulation compatible with $\varepsilon$ asymptotics by (1.2).

Let us stress that the extra cost of considering an augmented formulation is considered as affordable here because it is to be incorporated in a particle method. It is far from obvious that a similar strategy could be incorporated in different numerical methods.

The augmented system we choose here is

$$ (1.3) \qquad \left\{ \begin{aligned} \varepsilon \frac{d\mathbf{x}}{dt} &= \mathbf{w} - \varepsilon \left(e - \frac{1}{2}\|\mathbf{w}\|^2\right) \nabla_{\mathbf{x}}^{\perp} \left(\frac{1}{b}\right)(\mathbf{x}), \\ \varepsilon \frac{d\mathbf{e}}{dt} &= -\varepsilon \frac{d\phi(\mathbf{x})}{dt}, \\ \varepsilon \frac{d\mathbf{w}}{dt} &= \mathbf{E}(\mathbf{x}) - b(\mathbf{x}) \frac{\mathbf{w}^{\perp}}{\varepsilon}. \end{aligned} \right. $$

Now consider a time step $\Delta t > 0$ and define $(\mathbf{x}^n, e^n, \mathbf{w}^n)$, to serve as an approximation of $(\mathbf{x}, e, \mathbf{w})(n \Delta t)$ when $(\mathbf{x}, e, \mathbf{w})$ solves (1.3), through the scheme

$$ (1.4) \qquad \left\{ \begin{aligned} \frac{\mathbf{x}^{n+1} - \mathbf{x}^n}{\Delta t} &= \frac{\mathbf{w}^{n+1/2}}{\varepsilon} - \left( e^{n+1/2} - \frac{1}{2} \|\mathbf{w}^{n+1/2}\|^2 \right) \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\mathbf{x}^{n+1/2}), \\ \frac{e^{n+1} - e^n}{\Delta t} &= \frac{\phi(\mathbf{x}^n) - \phi(\mathbf{x}^{n+1})}{\Delta t}, \\ \frac{\mathbf{w}^{n+1} - \mathbf{w}^n}{\Delta t} &= \frac{1}{\varepsilon} \mathbf{E}(\mathbf{x}^{n+1/2}) - b(\mathbf{x}^{n+1/2}) \frac{(\mathbf{w}^{n+1/2})^{\perp}}{\varepsilon^2}, \end{aligned} \right. $$
where the underlining at half times denotes averaging

$$ \underline{\mathbf{z}}^{n+1/2} := \frac{\underline{\mathbf{z}}^{n+1} + \underline{\mathbf{z}}^n}{2}, $$

and initialization is $\mathbf{x}^0 = \mathbf{x}^0$, $\mathbf{w}^0 = \mathbf{v}^0$, $e^0 = \frac{1}{2}\|\mathbf{w}^0\|^2$. Note that by design (1.4) does not propagate the relation $e^n = \frac{1}{2}\|\mathbf{w}^n\|$ that holds initially. In Appendix A we check that Scheme (1.4) possesses a unique solution when $\Delta t$ is sufficiently small, independently of $\varepsilon$.

The scheme we consider here differs from the one studied in [9] mainly on the choice of the augmented system. Explicitly we have decided to rely on energy conservation to describe the kinetic energy evolution and to introduce in the evolution equation of $\mathbf{x}$ the extra "ghost" term — here $-\varepsilon(e - \frac{1}{2}\|\mathbf{w}\|^2) \nabla_{\mathbf{x}}^{\perp} (b^{-1})(\mathbf{x})$ — that vanishes identically at the continuous level for the kind of initial data we are interested in but contributes to enforce the right $\varepsilon$ asymptotics at the discrete level. In contrast, in the augmented system of [9] the dynamics of $(\mathbf{x}, e)$ is chosen to be

$$ \varepsilon \frac{d\mathbf{x}}{dt} = \mathbf{w}, $$

$$ \varepsilon \frac{de}{dt} = \mathbf{E}(\mathbf{x}) \cdot \mathbf{w}, $$

and a "ghost" term is introduced in the evolution equation of $\mathbf{w}$. Our numerical tests show similar performances between the scheme of [9] and the one studied here. The complete convergence study of the scheme of [9] would require a similar but lengthier analysis.

If we were to couple the underlying Vlasov equation with field equations we would need to reconstruct the velocity $\mathbf{v}^n$ from $(e^n, \mathbf{w}^n)$. A natural choice, compatible with energy conservation, is

$$ \mathbf{v}^n = \sqrt{2} e^n \frac{\mathbf{w}^n}{\|\mathbf{w}^n\|}, $$

when $e^n$ is nonnegative and $\mathbf{w}^n$ is non zero.

Even at the formal level it is not completely obvious that System (1.2) does describe the $\varepsilon$ asymptotics of the slow part of (1.1). In contrast at fixed $\Delta t$ one may easily guess what could be the asymptotic scheme for the slow part of solutions to (1.4) that do not blow up when $\varepsilon \to 0$. Indeed one may expect that at fixed $\Delta t$, for such solutions, $\underline{\mathbf{w}}^{n+1/2}\varepsilon^{-1} + (\underline{\mathbf{E}}^{\perp}/b)(\underline{\mathbf{x}}^{n+1/2})$ becomes $\varepsilon$-small as $\varepsilon$ goes to zero so that $(\mathbf{x}^n, e^n)$ becomes $\varepsilon$-close to $(\mathbf{y}^n, g^n)$ solving

$$ (1.5) \qquad \left\{ \begin{array}{l} \frac{\mathbf{y}^{n+1} - \mathbf{y}^n}{\Delta t} = -\frac{\underline{\mathbf{E}}^{\perp}}{b}(\underline{\mathbf{y}}^{n+1/2}) - g^{n+1/2} \nabla_{\mathbf{x}}^{\perp} \left(\frac{1}{b}\right) (\underline{\mathbf{y}}^{n+1/2}), \\ \frac{g^{n+1} - g^n}{\Delta t} = -\frac{\phi(\underline{\mathbf{y}}^{n+1}) - \phi(\underline{\mathbf{y}}^n)}{\Delta t}. \end{array} \right. $$

Note that (1.4) and (1.5) are second-order schemes respectively for (1.3) and (1.2). Such simple first guarantees are already provided in [5, 6, 8, 9] and are typically named there asymptotic consistency as $\varepsilon \to 0$ for a fixed $\Delta t$. Our present goal is to prove that indeed solutions of (1.4) do not blow up — inherently a stability question — and to provide uniform estimates when $(\varepsilon, \Delta t)$ is small but no relation between $\Delta t$ and (any power of) $\varepsilon$ is assumed.

**Main results.** Our main result is the following theorem that quantifies in which sense $(\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2)$ is indeed well-captured even when the time step $\Delta t$ is too large for oscillations at frequency $\varepsilon^{-2}$.

**Theorem 1.1.** Let $T > 0$, $M > 0$ and $\varepsilon_0 > 0$. There exist $\delta_0 > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, \mathbf{v})$ solves (1.1) on $[0, T]$ with initial datum $(\mathbf{x}^0, \mathbf{v}^0)$ such that $\|\mathbf{v}^0\| \le M$, and $(\mathbf{x}, e, \mathbf{w})$ solves (1.4) with initial datum $(\mathbf{x}^0, e^0, \mathbf{w}^0)$ such that $\mathbf{w}^0 = \mathbf{v}^0$, $e^0 = \frac{1}{2}\|\mathbf{v}^0\|^2$, then for any $n$ such that $0 \le n\Delta t \le T$,

$$ \|(\mathbf{x}^n, e^n) - (\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2) (n \Delta t)\| + \varepsilon \|\mathbf{w}^n - \mathbf{v}(n \Delta t)\| \le C \min \left( \left\{ \varepsilon + (\Delta t)^2 ; \frac{(\Delta t)^2}{\varepsilon^5} \right\} \right). $$
To give a more concrete outcome we point out that

$$ \min \left( \left\{ \varepsilon + (\Delta t)^2 ; \frac{(\Delta t)^2}{\varepsilon^5} \right\} \right) \le \begin{cases} \frac{(\Delta t)^2}{\varepsilon^5}, & \text{if } \Delta t \le \varepsilon^3, \quad \text{i.e. } (\Delta t)^{\frac{1}{3}} \le \varepsilon, \\ 2\varepsilon, & \text{if } \varepsilon^3 \le \Delta t \le \varepsilon^{\frac{1}{2}}, \quad \text{i.e. } (\Delta t)^2 \le \varepsilon \le (\Delta t)^{\frac{1}{3}}, \\ 2(\Delta t)^2, & \text{if } \varepsilon^{\frac{1}{2}} \le \Delta t, \quad \text{i.e. } \varepsilon \le (\Delta t)^2, \end{cases} $$

whereas

$$ \max_{0 < \varepsilon \le \varepsilon_0} \min \left( \left\{ \varepsilon + (\Delta t)^2 ; \frac{(\Delta t)^2}{\varepsilon^5} \right\} \right) \lesssim (\Delta t)^{\frac{1}{3}}. $$

Note that though the present analysis is more involved than the one of [10] the outcome is similar to the one obtained there for second-order schemes and, in particular, matches the best expectations rationalized in [10, Appendix A] for asymptotic preserving schemes of general oscillatory equations. Incidentally we stress that there are two main related ways in which our present analysis differ from the one in [10], both being responsible for new kinds of discretization consistency errors

* the inhomogeneity of the magnetic field forces more chain rule applications in the normal form part of the analysis, this being immaterial at the continuous level but resulting in consistency errors that are due to failure of exact discrete chain rules and which are denoted with $\tau_\#^\epsilon$ letters in discrete elimination lemmas, Lemmas 4.1, 4.2, 5.1 and 5.2

* estimates on $\varepsilon - \frac{1}{2}\|\mathbf{w}\|^2$ (see Lemma 5.5) should be thought as discretization consistency bounds and encode that the augmented formulation lets the scheme ensure the smallness of $\varepsilon - \frac{1}{2}\|\mathbf{w}\|^2$ when $\Delta t$ is smaller than $\varepsilon^2$ and its convergence to $\mathbf{g}$ otherwise.

The $\varepsilon + (\Delta t)^2$ part of the bound in Theorem 1.1 arises from

* the $\varepsilon$ proximity of $(\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2)(n\Delta t)$ and $(\mathbf{y}, g)(n\Delta t)$ (stated in Proposition 2.3)

* its discrete counterpart, the $\varepsilon$ proximity of $(\mathbf{x}^n, e^n)$ and $(\mathbf{y}^n, g^n)$ (stated in Proposition 4.3)

* the $(\Delta t)^2$ proximity of $(\mathbf{y}, g)(n\Delta t)$ and $(\mathbf{y}^n, g^n)$ which results from the fact that (1.5) is a second-order scheme for (1.2).

In turn the $(\Delta t)^2 \varepsilon^{-5}$ part is a "direct" convergence analysis encoding how the stiffness impacts the second-order convergence of (1.4) to (1.3) and the smallness of $\varepsilon - \frac{1}{2}\|\mathbf{w}\|^2$. Let us stress however that even the direct convergence analysis requires a fine understanding of the $\varepsilon$ asymptotics so as to tame the stiffness impact in the stability part of the convergence analysis. A rough stability analysis would provide a bound exhibiting a dramatic exponential growth as a power of $1/\varepsilon^2$.

In particular, even in the direct convergence part, our stability analysis involves the famous guiding-center variables $(\mathbf{x}_{gc}, e_{gc})$,

$$ \left\{ \begin{array}{l} \mathbf{x}_{gc} := \mathbf{x} - \varepsilon \frac{\mathbf{w}^{\perp}}{b(\mathbf{x})}, \\ e_{gc} := e + \frac{\varepsilon}{b(\mathbf{x})} \mathbf{E}^{\perp}(\mathbf{x}) \cdot \mathbf{w}. \end{array} \right. $$

These variables are $\varepsilon$ close to $(\mathbf{x}, e)$ but obey a slower dynamics. Indeed they remain $\varepsilon^2$ close to a suitable solution of (1.2); see Proposition 2.5. To keep our analysis as short as possible we have decided not include a convergence analysis of those. But we refer to [10] for an example of the latter and we point out that our numerical tests also show a better capture of $(\mathbf{x}_{gc}, e_{gc})$ following the rate

$$ \min \left( \left\{ \varepsilon^2 + (\Delta t)^2 ; \frac{(\Delta t)^2}{\varepsilon^4} \right\} \right) \le \begin{cases} \frac{(\Delta t)^2}{\varepsilon^4} & \text{if } \Delta t \le \varepsilon^3, \quad \text{i.e. } (\Delta t)^{\frac{1}{3}} \le \varepsilon, \\ 2\varepsilon^2 & \text{if } \varepsilon^3 \le \Delta t \le \varepsilon, \quad \text{i.e. } \Delta t \le \varepsilon \le (\Delta t)^{\frac{1}{3}}, \\ 2(\Delta t)^2 & \text{if } \varepsilon \le \Delta t, \end{cases} $$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Convergence errors, that is, comparison of the solution to the scheme (1.4) with the solution to (1.1). (Left) Comparison of the original slow variables $(\mathbf{x}, e)$. (Right) Comparison of their guiding center modifications $(\mathbf{x}_{gc}, e_{gc})$.</figcaption></figure>

which is consistent with the general discussion in [10, Appendix A]. Note that

$$ \max_{0 < \varepsilon \le \varepsilon_0} \min \left( \left\{ \varepsilon^2 + (\Delta t)^2 ; \frac{(\Delta t)^2}{\varepsilon^4} \right\} \right) \lesssim (\Delta t)^{\frac{2}{3}}. $$

**Numerical tests.** We conclude this introduction with numerical experiments confirming the sharpness of our estimates. In our experiments we choose

$$ b(\mathbf{x}) = \frac{10}{\sqrt{100 - \|\mathbf{x}\|^2}}, \qquad \phi(\mathbf{x}) = \frac{1}{2} \|\mathbf{x}\|^2. $$

Note that neither $b$ nor $\phi$ meet the assumptions of our results since $b$ is unbounded and not lower bounded away from zero and $\phi$ is not bounded with bounded derivatives. Their treatments would require slightly more elaborate arguments. Yet the outcome would be similar, as confirmed by numerical tests.

In the present experiments we fix the final time $T = 1$, the initial position to be $\mathbf{x}^0 = (2, 2)$ and the initial velocity to be $\mathbf{v}^0 = (3, 3)$. Moreover we constrain $\Delta t$ to be such that $N_{\Delta t} := T/(\Delta t)$ is an integer and measure sizes of vectors $(a_n)_{1 \le n \le N_{\Delta t}}$ through the discrete $\ell^1$ norm

$$ \|(a_n)_{1 \le n \le N_{\Delta t}}\|_{\Delta t} := \frac{1}{N_{\Delta t}} \sum_{n=1}^{N_{\Delta t}} \|a_n\|. $$

Numerical tests match quite closely theoretical bounds. For the sake of comparison, let us recall those here

* Figure 1, (Left) $\min(\{\varepsilon + (\Delta t)^2; (\Delta t)^2 \varepsilon^{-5}\})$, (Right) $\min(\{\varepsilon^2 + (\Delta t)^2; (\Delta t)^2 \varepsilon^{-4}\})$,

* Figure 2, (Left) $\varepsilon$, (Right) $\varepsilon^2$,

* Figure 3, (Left) $\varepsilon + (\Delta t)^2$, (Right) $\varepsilon^2 + (\Delta t)^2$.
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Discrete $\varepsilon$ asymptotic errors, that is, comparison of the solution to the scheme (1.4) with solutions to (1.5). (Left) Comparison of the original slow variables $(\mathbf{x}, e)$. (Right) Comparison of their guiding center modifications $(\mathbf{x}_{gc}, e_{gc})$. Note that initial data for (1.5) are set respectively to $(\mathbf{x}^0, e^0)$ for the computation of the left panel and to $(\mathbf{x}_{gc}^0, e_{gc}^0)$ for the computation of the right one.</figcaption></figure>

There are two departures worth commenting. On one hand, some of the regimes of the bounds reachable only with very small $\varepsilon$ are unseen. For instance this is the case of the $(\Delta t)^2$ regime for the bound on the convergence error in $(\mathbf{x}, e)$, which would be apparent if we were reaching the $\varepsilon \lesssim (\Delta t)^2$ regime. This numerical limitation of our experiments is due to the fact that computing reference solutions with very small $\varepsilon$ comes at a prohibitive price. On the other hand, for the convergence error in $(\mathbf{x}_{gc}, eGC)$, we measure $(\Delta t)^2 \varepsilon^{-3.3}$ in the regime where we should see $(\Delta t)^2 \varepsilon^{-4}$. This seems to be a purely numerical artifact.

**Outline.** After the present introduction the rest of the paper is organized as follows. In the next section we prove $\varepsilon$ asymptotics at the continuous level. After that, we begin the discrete analysis and provide at the discrete level, in this order, uniform bounds, $\varepsilon$ asymptotics and discrete error bounds. Then we combine those to prove Theorem 1.1. At last, in an appendix we prove that the implicit scheme may be solved without stiff constraint on the discretization.

*Acknowledgment*. KHT expresses his appreciation of the hospitality of IMT, Université Toulouse III, during the preparation of the present contribution.

## 2. CONTINUOUS ANALYSIS

Our continuous analysis follows closely the strategy of [7]: first prove uniform bounds then investigate the $\varepsilon$ asymptotics by tracking how multilinear expressions in the fast variable may be split into a slow part and a smaller fast part. However, in contrast with the analyses of [7], the restriction made on the structure provides readily relevant uniform bounds at the continuous level.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Continuous ε asymptotic errors, that is, comparison of the solution to the scheme (1.4) with solutions to (1.2). (Left) Comparison of the original slow variables (x, e). (Right) Comparison of their guiding center modifications (x<sub>gc</sub>, e<sub>gc</sub>). Note that initial data for (1.5) are set respectively to (x<sup>0</sup>, e<sup>0</sup>) for the computation of the left panel and to (x<sub>gc</sub><sup>0</sup>, e<sub>gc</sub><sup>0</sup>) for the computation of the right one.</figcaption></figure>

Indeed from the conservation of energy stems that solutions to (1.1) satisfy for any t and any ε > 0

$$
\|\mathbf{v}(t)\| \le \sqrt{\|\mathbf{v}^0\|^2 + 4 \|\phi\|_{L^\infty}} \le \|\mathbf{v}^0\| + 2 \sqrt{\|\phi\|_{L^\infty}}.
$$

2.1. **Original variables.** With the uniform bounds in hands we may turn to the elimination
process, aiming at a normal form where slow and fast variables are sufficiently uncoupled. The
keystone of the uncoupling is

$$
(2.1) \qquad \mathbf{v} = \epsilon^2 \frac{d}{dt} \left( \frac{\mathbf{v}^{\perp}}{b(\mathbf{x})} \right) - \epsilon \, d_{\mathbf{x}} \left( \frac{1}{b} \right) (\mathbf{x})(\mathbf{v}) \mathbf{v}^{\perp} - \epsilon \frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}),
$$

which is readily derived from (1.1). To begin with, inserting (2.1) back into the first equation of
(1.1) yields

$$
(2.2) \qquad \frac{d}{dt} \left( \mathbf{x} - \varepsilon \frac{\mathbf{v}^{\perp}}{b(\mathbf{x})} \right) = - d_{\mathbf{x}} \left( \frac{1}{b} \right) (\mathbf{x})(\mathbf{v}) \mathbf{v}^{\perp} - \frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}).
$$

At this stage, in order to proceed with the uncoupling, we need to eliminate quadratic terms.
Before doing so it is instructive to slightly generalize the elimination of **v** encoded by (2.1) as
elimination of linear terms. Combining (2.1) with direct differentiation proves the following lemma.

**Lemma 2.1.** *Consider L a differentiable function of time with values in linear maps from R² to some Rⁿ, n ∈ N*. Then solutions to (1.1) satisfy

$$
\mathbf{L}(t)(\mathbf{v}(t)) = -\epsilon^2 \frac{d}{dt} \frac{\kappa_{\mathbf{L}}(\cdot, \mathbf{x}, \mathbf{v})}{b(t)}(t) + \epsilon \eta_{\mathbf{L}}(t, \mathbf{x}(t), \mathbf{v}(t)),
$$
where

$$ \kappa_{\mathbf{L}}(t, \mathbf{x}, \mathbf{v}) := -\frac{1}{b(\mathbf{x})} \mathbf{L}(t)(\mathbf{v}^{\perp}), $$

and

$$ \eta_{\mathbf{L}}(t, \mathbf{x}, \mathbf{v}) := -d_{\mathbf{x}}\left(\frac{1}{b}\right)(\mathbf{x})(\mathbf{v})\mathbf{L}(t)(\mathbf{v}^{\perp}) - \frac{1}{b(\mathbf{x})}\mathbf{L}(t)(\mathbf{E}^{\perp}(\mathbf{x})) - \frac{\varepsilon}{b(\mathbf{x})}\frac{d}{dt}\mathbf{L}(t)(\mathbf{v}^{\perp}). $$

Above, to derive (2.2), we have just used the case when $\mathbf{L}$ is constant equal to the identity map.

Quadratic terms are not eliminated in the same way as linear terms are. They leave slow contributions. The starting observation in this direction is that if $\mathbf{B}$ is a bilinear form on $\mathbf{R}^2$ then for any vector $\mathbf{w}$

$$ (2.3) \qquad \mathbf{B}(\mathbf{w}, \mathbf{w}) + \mathbf{B}(\mathbf{w}^{\perp}, \mathbf{w}^{\perp}) = \|\mathbf{w}\|^2 \mathrm{Tr}(\mathbf{B}), $$

where Tr denotes the trace operator. Note also that for solutions to (1.1), $e := \frac{1}{2}\|\mathbf{v}\|^2$ is exactly as slow as $\mathbf{x}$ is, as may deduced from conservation of energy.

The following lemma provides a quadratic counterpart to Lemma 2.1

**Lemma 2.2.** Consider $\mathbf{B}$ a differentiable function of time with values in bilinear maps from $\mathbf{R}^2 \times \mathbf{R}^2$ to some $\mathbf{R}^n$, $n \in \mathbf{N}^*$. Then solutions to (1.1) satisfy

$$ \mathbf{B}(t)(\mathbf{v}(t), \mathbf{v}(t)) = e(t) \mathrm{Tr}(\mathbf{B}(t)) - \varepsilon^2 \frac{d}{dt} \kappa_{\mathbf{B}}(\cdot, \mathbf{x}, \mathbf{v}) (t) + \varepsilon \eta_{\mathbf{B}}(t, \mathbf{x}(t), \mathbf{v}(t)), $$

where

$$ \kappa_{\mathbf{B}}(t, \mathbf{x}, \mathbf{v}) := -\frac{1}{2} \frac{1}{b(\mathbf{x})} \mathbf{B}(t)(\mathbf{v}, \mathbf{v}^{\perp}), $$

and

$$ \eta_{\mathbf{B}}(t, \mathbf{x}, \mathbf{v}) := -\frac{1}{2} d_{\mathbf{x}}\left(\frac{1}{b}\right)(\mathbf{x})(\mathbf{v})\mathbf{B}(t)(\mathbf{v}, \mathbf{v}^{\perp}) - \frac{1}{2} \frac{\varepsilon}{b(\mathbf{x})} \frac{d}{dt}\mathbf{B}(t)(\mathbf{v}, \mathbf{v}^{\perp}) \\ - \frac{1}{2} \frac{1}{b(\mathbf{x})} \left(\mathbf{B}(t)(\mathbf{v}, \mathbf{E}^{\perp}(\mathbf{x})) + \mathbf{B}(t)(\mathbf{E}(\mathbf{x}), \mathbf{v}^{\perp})\right). $$

*Proof.* A direct computation relying again on (2.1) shows

$$ \frac{1}{2}\mathbf{B}(t)(\mathbf{v}(t), \mathbf{v}(t)) - \frac{1}{2}\mathbf{B}(t)(\mathbf{v}(t)^{\perp}, \mathbf{v}(t)^{\perp}) = -\varepsilon^2 \frac{d}{dt} \kappa_{\mathbf{B}}(\cdot, \mathbf{x}, \mathbf{v}) (t) + \varepsilon \eta_{\mathbf{B}}(t, \mathbf{x}(t), \mathbf{v}(t)), $$

whereas (2.3) implies

$$ \mathbf{B}(t)(\mathbf{v}(t), \mathbf{v}(t)) = e(t) \mathrm{Tr}(\mathbf{B}(t)) + \frac{1}{2}\mathbf{B}(t)(\mathbf{v}(t), \mathbf{v}(t)) - \frac{1}{2}\mathbf{B}(t)(\mathbf{v}(t)^{\perp}, \mathbf{v}(t)^{\perp}). $$

□

Coming back to (2.2) we apply Lemma 2.2 with $\mathbf{B}_1(t) := \mathcal{B}_1(\mathbf{x}(t))$

$$ \mathcal{B}_1(\mathbf{x})(\mathbf{w}_1, \mathbf{w}_2) := -d_{\mathbf{x}}\left(\frac{1}{b}\right)(\mathbf{x})(\mathbf{w}_1)\mathbf{w}_2^{\perp} $$

and observe that for this bilinear form setting $\mathbf{w} = \nabla_{\mathbf{x}}(b^{-1})(\mathbf{x}(t))$ in (2.3) yields

$$ \mathrm{Tr}(\mathbf{B}_1(t)) = -\nabla_{\mathbf{x}}^{\perp}\left(\frac{1}{b}\right)(\mathbf{x}(t)). $$
Therefore solutions to (1.1) satisfy

$$
(2.4) \qquad \frac{d}{dt} (\mathbf{x} + \varepsilon \kappa_1(\mathbf{x}, \mathbf{v})) = -e \nabla_{\mathbf{x}}^{\perp} \left(\frac{1}{b}\right) (\mathbf{x}) - \frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}) + \varepsilon \eta_1(\mathbf{x}, \mathbf{v}),
$$

with

$$
\kappa_1(\mathbf{x}, \mathbf{v}) := - \frac{\mathbf{v}^{\perp}}{b(\mathbf{x})} - \frac{\varepsilon}{2} \frac{1}{b(\mathbf{x})} \mathcal{B}_1(\mathbf{x})(\mathbf{v}, \mathbf{v}^{\perp}),
$$

and

$$
\eta_1(\mathbf{x}, \mathbf{v}) := \frac{1}{4} d_{\mathbf{x}}^2 \left( \frac{1}{b^2} \right) (\mathbf{x})(\mathbf{v}, \mathbf{v}) \mathbf{v} - \frac{1}{2} \frac{1}{b(\mathbf{x})} \left( \mathcal{B}_1(\mathbf{x})(\mathbf{v}, \mathbf{E}^{\perp}(\mathbf{x})) + \mathcal{B}_1(\mathbf{x})(\mathbf{E}(\mathbf{x}), \mathbf{v}^{\perp}) \right).
$$

We have reached the stage when we may derive the ε asymptotic behavior of (x, e).

**Proposition 2.3.** Let $T > 0$, $M > 0$ and $\varepsilon_0 > 0$. There exists $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ if $(\mathbf{x}, \mathbf{v})$ solves (1.1) with initial datum $(\mathbf{x}^0, \mathbf{v}^0)$ such that $\|\mathbf{v}^0\| \le M$ and $(\mathbf{y}, g)$ solves (1.2) with initial datum $(\mathbf{y}^0, g^0) = (\mathbf{x}^0, \frac{1}{2}\|\mathbf{v}^0\|^2)$ then with $e := \frac{1}{2}\|\mathbf{v}\|^2$ for any $0 \le t \le T$,

$$
\|(\mathbf{x}(t), e(t)) - (\mathbf{y}(t), g(t))\| \le C\varepsilon.
$$

The dependence of C on T and ||v^0|| could actually be made rather explicit ; see the related
analyses in [7, 10] but we have decided to cut technicalities to a bare minimum.

Proof. Note that from the conservation of energy one deduces that for some $M'$ depending only on $M$ (and $\phi$), $e$, $\mathbf{v}$ and $g$ remain bounded by $M'$ when $\|\mathbf{v}^0\| \le M$. Then observe that

$$
(\mathbf{z}, \mathbf{w}) \mapsto \kappa_1(\mathbf{z}, \mathbf{w}), \qquad (\mathbf{z}, \mathbf{w}) \mapsto \eta_1(\mathbf{z}, \mathbf{w})
$$

are uniformly bounded and

$$
(\mathbf{z}, h) \mapsto -h \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\mathbf{z}) - \frac{\mathbf{E}^{\perp}}{b} (\mathbf{z})
$$

is uniformly Lipschitz when restricted respectively to $\|\mathbf{w}\| \le M'$ and $|h| \le M'$.

By integrating (2.4) one deduces that for some $K$ depending only on $M > 0$ and $\varepsilon_0 > 0$ (and the fields), for any $0 \le s \le t$,

$$
\| \mathbf{x}(t) - \mathbf{x}(s) + \int_s^t \left( e(\tau) \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\mathbf{x}(\tau)) + \frac{\mathbf{E}^{\perp}}{b} (\mathbf{x}(\tau)) \right) d\tau \| \le K \varepsilon (1+t-s).
$$

In turn note that for any $0 \le s \le t$,

$$
\mathbf{y}(t) - \mathbf{y}(s) + \int_s^t \left( g(\tau) \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\mathbf{y}(\tau)) + \frac{\mathbf{E}^{\perp}}{b} (\mathbf{y}(\tau)) \right) d\tau
$$

and for any $t$, $e(t) - g(t) = -(\phi(\mathbf{x}(t)) - \phi(\mathbf{y}(t)))$. Subtracting and using the aforementioned Lipschitz bound yield that for some $K'$ depending only on $M > 0$ and $\varepsilon_0 > 0$ (and the fields), for any $t \ge 0$,

$$
\| \mathbf{y}(t) - \mathbf{x}(t) \| \le K' \int_0^t \| \mathbf{y}(s) - \mathbf{x}(s) \| \, d\,s + K \varepsilon (1+t).
$$

The proof is then concluded by the Grönwall lemma with $C = e^{K'T} K (1+T) (1 + \|\nabla_x \phi\|_{L^\infty})$ since
the bound on $e-g$ stems from the bound on $\mathbf{x}-\mathbf{y}$ and the energy conservation. $\square$
## 2.2. Guiding-center variables.
The above computations also suggest that the evolution of the guiding center variable

$$x_{gc} := x - \varepsilon \frac{v^{\perp}}{b(x)}$$

could be described as part of a slow uncoupled dynamics up to $O(\varepsilon^2)$ errors.

In order to achieve this we need to refine (2.4) and the energy conservation so as to get rid of $O(\varepsilon)$ terms involving $v$ either appearing explicitly in $\varepsilon \eta_1$ or arising from the replacement of $x$ with $x_{gc}$ in $O(1)$ terms. Heeding the shape of $\eta_1$ we need a trilinear version of Lemmas 2.1 and 2.2.

**Lemma 2.4.** Consider $\mathbf{T}$ a differentiable function of time with values in symmetric trilinear maps from $\mathbf{R}^2 \times \mathbf{R}^2 \times \mathbf{R}^2$ to some $\mathbf{R}^n$, $n \in \mathbf{N}^*$. Then solutions to (1.1) satisfy

$$\mathbf{T}(t)(\mathbf{v}(t), \mathbf{v}(t), \mathbf{v}(t)) = -\varepsilon^2 \frac{d \kappa_{\mathbf{T}}(\cdot, \mathbf{x}, \mathbf{v})}{dt}(t) + \varepsilon \eta_{\mathbf{T}}(t, \mathbf{x}(t), \mathbf{v}(t)),$$

where

$$\kappa_{\mathbf{T}}(t, \mathbf{x}, \mathbf{v}) := - \frac{1}{b(\mathbf{x})} \left( \mathbf{T}(t)(\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) + \frac{2}{3} \mathbf{T}(t)(\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp}) \right),$$

and

$$\eta_{\mathbf{T}}(t, \mathbf{x}, \mathbf{v}) := - d_{\mathbf{x}} \left(\frac{1}{b}(\mathbf{x})(\mathbf{v})\right) \left(\mathbf{T}(t)(\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) + \frac{2}{3} \mathbf{T}(t)(\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp})\right) \\ - \frac{\varepsilon}{b(\mathbf{x})} \left( \frac{d \mathbf{T}}{dt}(t)(\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) + \frac{2}{3} \frac{d \mathbf{T}}{dt}(t)(\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp}) \right) \\ - \frac{1}{b(\mathbf{x})} \left( 2 \mathbf{T}(t)(\mathbf{E}(\mathbf{x}), \mathbf{v}, \mathbf{v}^{\perp}) + \mathbf{T}(t)(\mathbf{E}^{\perp}(\mathbf{x}), \mathbf{v}, \mathbf{v}) + 2 \mathbf{T}(t)(\mathbf{E}^{\perp}(\mathbf{x}), \mathbf{v}^{\perp}, \mathbf{v}^{\perp}) \right).$$

Let us set

$$e_{gc} := e + \frac{\varepsilon}{b(\mathbf{x})} d_{\mathbf{x}} \phi(\mathbf{x}) (\mathbf{v}^{\perp}) = e + \frac{\varepsilon}{b(\mathbf{x})} \mathbf{E}^{\perp}(\mathbf{x}) \cdot \mathbf{v}$$

so that solutions to (1.1) satisfy

$$(2.5) \qquad \frac{d}{dt} (e_{gc} + \phi(\mathbf{x}_{gc}) + \varepsilon^2 \kappa_e(\mathbf{x}, \mathbf{v})) = 0$$

with

$$\kappa_e(\mathbf{x}, \mathbf{v}) := - \frac{1}{b(\mathbf{x})^2} \int_0^1 d_{\mathbf{x}}^2 \phi \left( \mathbf{x} - \varepsilon \tau \frac{\mathbf{v}^{\perp}}{b(\mathbf{x})} \right) (\mathbf{v}^{\perp}, \mathbf{v}^{\perp}) (1-\tau) d\tau.$$

Now we apply Lemma 2.1 with $\mathbf{L}_2(t) := \mathcal{L}_2(\mathbf{x}(t), \mathbf{x}_{gc}(t), e(t))$ where

$$\mathcal{L}_2(\mathbf{x}, \mathbf{y}, e)(\mathbf{w}) := - \frac{e}{b(\mathbf{x})} d \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \right) (\mathbf{x})(\mathbf{w}^{\perp}) + \frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}) \cdot \mathbf{w} \quad \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\mathbf{y}) \\ - \frac{1}{b(\mathbf{x})} d_{\mathbf{x}} \left( \frac{\mathbf{E}^{\perp}}{b} \right) (\mathbf{x})(\mathbf{w}^{\perp}) - \frac{1}{2} \frac{1}{b(\mathbf{x})} \left( \mathcal{B}_1(\mathbf{x})(\mathbf{w}, \mathbf{E}^{\perp}(\mathbf{x})) + \mathcal{B}_1(\mathbf{x})(\mathbf{E}(\mathbf{x}), \mathbf{w}^{\perp}) \right)$$

and Lemma 2.4 with $\mathbf{T}_1(t) := \mathcal{T}_1(\mathbf{x}(t))$ where

$$\mathcal{T}_1(\mathbf{x})(\mathbf{w}_1, \mathbf{w}_2, \mathbf{w}_3) \\ := \frac{1}{12} \left( d_{\mathbf{x}}^2 \left( \frac{1}{b^2} \right) (\mathbf{x})(\mathbf{w}_1, \mathbf{w}_2) \mathbf{w}_3 + d_{\mathbf{x}}^2 \left( \frac{1}{b^2} \right) (\mathbf{x})(\mathbf{w}_3, \mathbf{w}_1) \mathbf{w}_2 + d_{\mathbf{x}}^2 \left( \frac{1}{b^2} \right) (\mathbf{x})(\mathbf{w}_2, \mathbf{w}_3) \mathbf{w}_1 \right).$$
This implies that solutions to (1.1) satisfy

$$ (2.6) \quad \frac{d}{dt} (\mathbf{x}_{gc} + \varepsilon^2 \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{v})) = -e_{gc} \nabla_{\mathbf{x}}^{\perp} \left(\frac{1}{b}\right) (\mathbf{x}_{gc}) - \frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}_{gc}) + \varepsilon^2 \eta_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{v}), $$

with $\kappa_2$ given by

$$ \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{v}) := -\frac{1}{2} \frac{1}{b(\mathbf{x})} \mathcal{B}_1(\mathbf{x})(\mathbf{v}, \mathbf{v}^{\perp}) \\ - \frac{\varepsilon}{b(\mathbf{x})} (\mathcal{L}_2(\mathbf{x}, \mathbf{x}_{gc}, e)(\mathbf{v}^{\perp}) + \mathcal{T}_1(\mathbf{x})(\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) + \frac{2}{3} \mathcal{T}_1(\mathbf{x})(\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp})). $$

whereas $\eta_2$ is

$$ \eta_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{v}) := \frac{1}{b(\mathbf{x})^2} \int_0^1 d_{\mathbf{x}}^2 \left( e \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) + \frac{\mathbf{E}^{\perp}}{b} \right) \left( \mathbf{x} - \varepsilon \tau \frac{\mathbf{v}^{\perp}}{b(\mathbf{x})} \right) (\mathbf{v}^{\perp}, \mathbf{v}^{\perp}) (1-\tau) d\tau \\ - \left( d_{\mathbf{x}} \left( \frac{\mathcal{L}_2}{b} \right) (\mathbf{x}, \mathbf{x}_{gc}, e)(\mathbf{v}) \right) (\mathbf{v}^{\perp}) - \frac{\mathcal{L}_2(\mathbf{x}, \mathbf{x}_{gc}, e)}{b(\mathbf{x})} (\mathbf{E}^{\perp}(\mathbf{x})) \\ - \frac{1}{b(\mathbf{x})} \left( d_e \mathcal{L}_2(\mathbf{x}, \mathbf{x}_{gc}, e) (\mathbf{E}(\mathbf{x}) \cdot \mathbf{v}) \right) (\mathbf{v}^{\perp}) \\ + \frac{\varepsilon}{b(\mathbf{x})} \left( d_y \mathcal{L}_2(\mathbf{x}, \mathbf{x}_{gc}, e) \left( d_{\mathbf{x}} \left( \frac{1}{b} \right) (\mathbf{x})(\mathbf{v}) \mathbf{v}^{\perp} + \frac{\mathbf{E}^{\perp}}{b}(\mathbf{x}) \right) \right) (\mathbf{v}^{\perp}) \\ - \left( \left( d_{\mathbf{x}} \left( \frac{\mathcal{T}_1}{b} \right) (\mathbf{x})(\mathbf{v}) \right) (\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) + \frac{2}{3} \left( d_{\mathbf{x}} \left( \frac{\mathcal{T}_1}{b} \right) (\mathbf{x})(\mathbf{v}) \right) (\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp}) \right) \\ - \left( 2 \frac{\mathcal{T}_1}{b}(\mathbf{x})(\mathbf{E}(\mathbf{x}), \mathbf{v}, \mathbf{v}^{\perp}) + \frac{\mathcal{T}_1}{b}(\mathbf{x})(\mathbf{E}^{\perp}(\mathbf{x}), \mathbf{v}, \mathbf{v}) + 2 \frac{\mathcal{T}_1}{b}(\mathbf{x})(\mathbf{E}^{\perp}(\mathbf{x}), \mathbf{v}^{\perp}, \mathbf{v}^{\perp}) \right). $$

From here a slight variation on the proof of Proposition 2.3 proves the following.

**Proposition 2.5.** Let $T > 0$, $M > 0$ and $\varepsilon_0 > 0$. There exists $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ if $(\mathbf{x}, \mathbf{v})$ solves (1.1) with initial datum $(\mathbf{x}^0, \mathbf{v}^0)$ such that $\|\mathbf{v}^0\| \le M$ and $(\mathbf{y}, g)$ solves (1.2) with initial datum

$$ (\mathbf{y}^0, g^0) = \left( \mathbf{x}^0 - \frac{\varepsilon}{b(\mathbf{x}^0)} (\mathbf{v}^0)^{\perp}, \quad \frac{1}{2} \|\mathbf{v}^0\|^2 + \frac{\varepsilon}{b(\mathbf{x}^0)} \mathbf{E}^{\perp}(\mathbf{x}^0) \cdot \mathbf{v}^0 \right) $$

then with

$$ \mathbf{x}_{gc} := \mathbf{x} - \frac{\varepsilon}{b(\mathbf{x})} \mathbf{v}^{\perp}, \qquad e_{gc} := \frac{1}{2} \|\mathbf{v}\|^2 + \frac{\varepsilon}{b(\mathbf{x})} \mathbf{E}^{\perp}(\mathbf{x}) \cdot \mathbf{v}, $$

for any $0 \le t \le T$,

$$ \|(\mathbf{x}_{gc}(t), e_{gc}(t)) - (\mathbf{y}(t), g(t))\| \le C \varepsilon^2. $$

### 3. DISCRETE UNIFORM BOUNDS

From now on we shall focus on discrete analysis. It is useful to introduce a few more pieces of notation so as to gain in writing concision. We shall denote composition with a subscript

$$ \mathbf{G}_{\mathbf{z}}^{\ell} := \mathbf{G}(\mathbf{z}^{\ell}) $$

so that for instance

$$ \mathbf{G}_{\mathbf{z}}^{n+1/2} = \mathbf{G} \left( \frac{\mathbf{z}^{n+1} + \mathbf{z}^n}{2} \right), \qquad \underline{\mathbf{G}_{\mathbf{z}}}^{n+1/2} = \frac{\mathbf{G}(\mathbf{z}^{n+1}) + \mathbf{G}(\mathbf{z}^n)}{2}. $$
Moreover we introduce the matrix $\mathbf{J}$ associated with $\mathbf{J}^{\perp}$

$$ (3.1) \qquad \mathbf{J} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}, $$

so that $\mathbf{z}^{\perp} = \mathbf{J}\mathbf{z}$. $\mathbf{J}$ is skew-adjoint and unitary. In particular $\mathbf{J}^{-1} = -\mathbf{J}$. For concision's sake it is also useful to introduce

$$ \mathbf{F} := -\mathbf{b}^{-1}\mathbf{J}\mathbf{E}. $$

The following lemma gathers key observations on $\mathbf{J}$ towards uniform bounds and stability.

**Lemma 3.1.** (1) For any real $\alpha$, $\mathrm{Id} + \alpha \mathbf{J}$ is invertible and for any vector $\mathbf{z}$

$$ \|(\mathrm{Id} + \alpha \mathbf{J}) \mathbf{z}\| = \sqrt{1 + \alpha^2} \|\mathbf{z}\|, \qquad \|(\mathrm{Id} + \alpha \mathbf{J})^{-1} \mathbf{z}\| = \frac{1}{\sqrt{1 + \alpha^2}} \|\mathbf{z}\|. $$

(2) For any real $\alpha$, $(\mathrm{Id} + \alpha \mathbf{J})^{-1}(\mathrm{Id} - \alpha \mathbf{J})$ is unitary.

(3) For any real $\alpha, \alpha'$,

$$ (\mathrm{Id} + \alpha \mathbf{J})^{-1} - (\mathrm{Id} + \alpha' \mathbf{J})^{-1} = -(\mathrm{Id} + \alpha \mathbf{J})^{-1} (\alpha - \alpha') \mathbf{J} (\mathrm{Id} + \alpha' \mathbf{J})^{-1}. $$

*Proof.* The first norm equality stems from the Pythagoras' equality and the fact that $\mathbf{J}$ is skew-adjoint and unitary. It implies that $\mathrm{Id} + \alpha \mathbf{J}$ is one-to-one hence invertible, the second norm equality and the unitary property.

The comparison equality may be derived by composing from the left with $(\mathrm{Id} + \alpha \mathbf{J})^{-1}$ and from the right with $(\mathrm{Id} + \alpha' \mathbf{J})^{-1}$ the simple $(\mathrm{Id} + \alpha' \mathbf{J}) - (\mathrm{Id} + \alpha \mathbf{J}) = -(\alpha - \alpha')\mathbf{J}$. $\square$

The goal of the present section is to prove the following proposition.

**Proposition 3.2.** Let $T > 0$, $M > 0$, $\varepsilon_0 > 0$. There exist $\delta_0 > 0$, $M' > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, e, \mathbf{w})$ solves (1.4) with initial datum $(\mathbf{x}^0, e^0, \mathbf{w}^0)$ such that $\|\mathbf{w}^0\| \le M$, $|e^0| \le \frac{1}{2}M^2$ then when $0 \le n \Delta t \le T$,

$$ |e^n| + \|\mathbf{w}^n\| \le M', $$

and

$$ \left\{ \begin{array}{l} \|\underline{\mathbf{w}}^{n+1/2}\| \le C \left( \frac{1}{\sqrt{1+\lambda^2}} + \varepsilon \right), \\ \|\mathbf{x}^{n+1} - \mathbf{x}^n\| + \varepsilon \|\mathbf{w}^{n+1} - \mathbf{w}^n\| \le C \min \left( \left\{ \frac{\Delta t}{\varepsilon} ; \Delta t + \varepsilon \right\} \right), \end{array} \right. $$

where $\lambda := \Delta t / \varepsilon^2$.

A closer inspection of the proof, which we omit for the sake of readability, reveals that $\delta_0$ depends on the fields but not on $(T, M, \varepsilon_0)$.

Let us begin the proof by observing that as in the continuous case the bound on $e^n$ follows readily from the energy conservation provided that $M' \ge \frac{1}{2}M^2 + 2\|\phi\|_{\infty}$. The proofs of the rest of the bounds require more work and are scattered in the following subsections.

**Bound on $\underline{\mathbf{w}}^{n+1/2}$.** To begin with we set for any $n \ge 0$

$$ \mathbf{z}^n := \underline{\mathbf{w}}^{n+1/2} - \varepsilon \mathbf{F}_{\underline{\mathbf{x}}}^{n+1/2} $$

so that

$$ \mathbf{w}^{n+1} - \mathbf{w}^n = -\lambda b_{\underline{\mathbf{x}}}^{n+1/2} \mathbf{J} \mathbf{z}^n. $$

Since for any $n \ge 0$, $\mathbf{w}^{n+2} - \mathbf{w}^n = 2(\underline{\mathbf{w}}^{n+3/2} - \underline{\mathbf{w}}^{n+1/2})$, we deduce by subtracting two consecutive steps that for any $n \ge 0$

$$ \mathbf{z}^{n+1} - \mathbf{z}^n = -\varepsilon (\mathbf{F}_{\underline{\mathbf{x}}}^{n+3/2} - \mathbf{F}_{\underline{\mathbf{x}}}^{n+1/2}) - \frac{1}{2}\lambda b_{\underline{\mathbf{x}}}^{n+3/2} \mathbf{J} \mathbf{z}^{n+1} + \frac{1}{2}\lambda b_{\underline{\mathbf{x}}}^{n+1/2} \mathbf{J} \mathbf{z}^n, $$
or in other words

$$ (\mathrm{Id} + \frac{1}{2}\lambda b_{\underline{x}}^{n+3/2} \mathbf{J}) \mathbf{z}^{n+1} = (\mathrm{Id} + \frac{1}{2}\lambda b_{\underline{x}}^{n+1/2} \mathbf{J}) \mathbf{z}^{n} - \varepsilon (\mathbf{F}_{\underline{x}}^{n+3/2} - \mathbf{F}_{\underline{x}}^{n+1/2}). $$

At this stage, Lemma 3.1 motivates the introduction of

$$ \zeta^n := \sqrt{1 + \frac{1}{4}\lambda^2 (b_{\underline{x}}^{n+1/2})^2} \|\mathbf{z}^n\| $$

so that for any $n \ge 0$

$$ \zeta^{n+1} \le \zeta^n + \varepsilon \|\mathbf{F}_{\underline{x}}^{n+3/2} - \mathbf{F}_{\underline{x}}^{n+1/2}\|. $$

To proceed we observe that

$$ \begin{aligned} \varepsilon \|\mathbf{F}_{\underline{x}}^{n+3/2} - \mathbf{F}_{\underline{x}}^{n+1/2}\| & \lesssim \varepsilon \|\underline{x}^{n+3/2} - \underline{x}^{n+1/2}\| \lesssim \varepsilon \|\underline{x}^{n+2} - \underline{x}^{n+1}\| + \varepsilon \|\underline{x}^{n+1} - \underline{x}^n\| \\ & \lesssim \Delta t \left( \|\mathbf{z}^{n+1}\| + \|\mathbf{z}^n\| + \varepsilon + \varepsilon \|\underline{\mathbf{w}}^{n+3/2}\|^2 + \varepsilon \|\underline{\mathbf{w}}^{n+1/2}\|^2 \right) \\ & \lesssim \Delta t \left( \zeta^{n+1} + \zeta^n + \varepsilon + \varepsilon \left( \|\mathbf{w}^{n+2}\|^2 + \|\mathbf{w}^{n+1}\|^2 + \|\mathbf{w}^n\|^2 \right) \right). \end{aligned} $$

To complete the analysis we study now the variations of $\varepsilon \|\mathbf{w}^n\|^2$. To do so we first observe that by taking the scalar product of the $\mathbf{w}$ equation of (1.4) with $\underline{\mathbf{w}}^{n+1/2}$ one obtains

$$ \varepsilon \|\mathbf{w}^{n+1}\|^2 - \varepsilon \|\mathbf{w}^n\|^2 = 2 \Delta t \mathbf{E}_{\underline{x}}^{n+1/2} \cdot \underline{\mathbf{w}}^{n+1/2}, $$

so that

$$ \varepsilon \|\mathbf{w}^{n+1}\|^2 - \varepsilon \|\mathbf{w}^n\|^2 \lesssim \Delta t (\zeta^n + \varepsilon) $$

Introducing for any $n$

$$ Z^n := \zeta^n + \varepsilon + \varepsilon \|\mathbf{w}^{n+1}\|^2 + \varepsilon \|\mathbf{w}^n\|^2 $$

and gathering the above estimates we deduce that for some constant $K$ and any $n$

$$ Z^{n+1} \le Z^n + K \Delta t (Z^{n+1} + Z^n). $$

Assuming that $K\Delta t \le 1/2$, one derives that for any $n$

$$ Z^{n+1} \le \frac{1+K\Delta t}{1-K\Delta t} Z^n = \left(1+\frac{2K\Delta t}{1-K\Delta t}\right) Z^n \le (1+4K\Delta t) Z^n. $$

By iteration one deduces that for any $n$

$$ \|\underline{\mathbf{w}}^{n+1/2}\| \lesssim \|\underline{\mathbf{z}}^n\| + \varepsilon \lesssim \frac{Z^n}{\sqrt{1+\lambda^2}} + \varepsilon \lesssim \frac{(1+4K\Delta t)^n Z^0}{\sqrt{1+\lambda^2}} + \varepsilon \lesssim e^{4K n \Delta t} (\|\underline{\mathbf{w}}^{1/2}\| + \varepsilon). $$

To conclude the bound on $\|\underline{\mathbf{w}}^{n+1/2}\|$, there only remains to estimate $\|\underline{\mathbf{w}}^{1/2}\|$. Once again we go back to the $\mathbf{w}$ equation of (1.4) and deduce this time

$$ \underline{\mathbf{w}}^{1/2} = (\mathrm{Id} + \frac{1}{2}\lambda b_{\underline{x}}^{1/2} \mathbf{J})^{-1} \left( \mathbf{w}^0 + \frac{1}{2}\varepsilon\lambda \mathbf{E}_{\underline{x}}^{1/2} \right). $$

Thus, as expected,

$$ \|\underline{\mathbf{w}}^{1/2}\| \lesssim \frac{1}{\sqrt{1+\lambda^2}} + \varepsilon $$

and the claimed bound on $\|\underline{\mathbf{w}}^{n+1/2}\|$ holds.
**Bound on w**<sup>n</sup>. For the next step of the proof the starting point is the w equation of (1.4) written as

$$ (\mathrm{Id} - \frac{1}{2}\lambda b_{\underline{x}}^{n+1/2} \mathbf{J}) \mathbf{w}^{n+1} = (\mathrm{Id} + \frac{1}{2}\lambda b_{\underline{x}}^{n+1/2} \mathbf{J}) \mathbf{w}^{n} - \varepsilon \lambda b_{\underline{x}}^{n+1/2} \mathbf{J} \mathbf{F}_{\underline{x}}^{n+1/2}. $$

Introducing for any $n \ge 1$

$$ \mathbf{u}^n := \mathbf{w}^n - \varepsilon \mathbf{J} \mathbf{F}_{\underline{x}}^{n-1/2} $$

one deduces that for any $n \ge 1$

$$ (\mathrm{Id} - \frac{1}{2}\lambda b_{\underline{x}}^{n+1/2} \mathbf{J}) \mathbf{u}^{n+1} = (\mathrm{Id} + \frac{1}{2}\lambda b_{\underline{x}}^{n+1/2} \mathbf{J}) \mathbf{u}^{n} - \varepsilon (\mathrm{Id} + \frac{1}{2}\lambda b_{\underline{x}}^{n+1/2} \mathbf{J}) (\mathbf{F}_{\underline{x}}^{n+1/2} - \mathbf{F}_{\underline{x}}^{n-1/2}). $$

so that

$$ \|\mathbf{u}^{n+1}\| \le \|\mathbf{u}^n\| + \varepsilon \|\mathbf{F}_{\underline{x}}^{n+1/2} - \mathbf{F}_{\underline{x}}^{n-1/2}\|. $$

Benefiting from the bounds of the previous step we obtain for any $n \ge 1$ such that $n\Delta t \le T$

$$ \|\mathbf{u}^{n+1}\| - \|\mathbf{u}^n\| \lesssim \Delta t \left( \frac{1}{\sqrt{1 + \lambda^2}} + \varepsilon \right). $$

By iteration this gives for any $n \ge 1$ such that $n\Delta t \le T$

$$ \begin{aligned} \|\mathbf{w}^n\| & \lesssim \|\mathbf{u}^n\| + \varepsilon \lesssim \|\mathbf{u}^1\| + \varepsilon + n\Delta t \left( \frac{1}{\sqrt{1 + \lambda^2}} + \varepsilon \right) \\ & \lesssim \|\mathbf{w}^1\| + \frac{1}{\sqrt{1 + \lambda^2}} + \varepsilon \lesssim \|\mathbf{w}^0\| + \frac{1}{\sqrt{1 + \lambda^2}} + \varepsilon \lesssim 1, \end{aligned} $$

$$ \text{since } \mathbf{w}^1 = \mathbf{w}^0 + 2\underline{\mathbf{w}}^{1/2}. $$

**Bounds on x<sup>n+1</sup> - x<sup>n</sup> and w<sup>n+1</sup> - w<sup>n</sup>.** The last step follows readily from

$$ \begin{aligned} \mathbf{x}^{n+1} - \mathbf{x}^n &= \varepsilon \lambda \underline{\mathbf{w}}^{n+1/2} - \Delta t \left( \frac{\varepsilon^{n+1/2}}{2} - \frac{1}{2} \|\underline{\mathbf{w}}^{n+1/2}\|^2 \right) \nabla_{\underline{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{x}^{n+1/2}) \\ \varepsilon (\mathbf{w}^{n+1} - \mathbf{w}^n) &= -\varepsilon \lambda b_{\underline{x}}^{n+1/2} \mathbf{J} \underline{\mathbf{w}}^{n+1/2} + \Delta t \mathbf{E}_{\underline{x}}^{n+1/2} \end{aligned} $$

and the fact that when $0 \le n\Delta t \le T$

$$ \varepsilon \lambda \|\underline{\mathbf{w}}^{n+1/2}\| \lesssim \varepsilon \frac{\lambda}{\sqrt{1 + \lambda^2}} + \varepsilon^2 \lambda \lesssim \min \left( \left\{ \frac{\Delta t}{\varepsilon}; \varepsilon \right\} \right) + \Delta t. $$

## 4. ASYMPTOTIC ANALYSIS

The goal of the present section is to provide a discrete counterpart to Proposition 2.3. With uniform bounds in hands, the main ingredients to add are multilinear elimination lemmas.

To prepare those we first gather bilinear discrete differentiation and averaging rules

$$ \left\{ \begin{aligned} (ab)^{n+1} - (ab)^n &= \underline{a}^{n+1/2} (b^{n+1} - b^n) + (a^{n+1} - a^n) \underline{b}^{n+1/2}, \\ (\underline{ab})^{n+1/2} &= \underline{a}^{n+1/2} \underline{b}^{n+1/2} + \frac{1}{4} (a^{n+1} - a^n) (b^{n+1} - b^n). \end{aligned} \right. $$

From (1.4) stems the following discrete counterpart to (2.1)

$$ (4.1) \qquad \underline{\mathbf{w}}^{n+1/2} = \frac{\varepsilon^2}{\Delta t} \frac{(\mathbf{w}^{n+1})^{\perp} - (\mathbf{w}^n)^{\perp}}{b_{\underline{x}}^{n+1/2}} - \varepsilon \left( \frac{\mathbf{E}^{\perp}}{b} \right)_{\underline{x}}^{n+1/2}. $$

In turn, straightforward computations imply the following discrete counterparts to Lemmas 2.1 and 2.2.
**Lemma 4.1.** Consider $\mathcal{L}$ a differentiable function of space with values in linear maps from $\mathbb{R}^2$ to some $\mathbb{R}^n$, $n \in \mathbb{N}^*$. Then solutions to (1.4) satisfy

$$ (\mathcal{L}_{\underline{\mathbf{x}}}(\underline{\mathbf{w}}))^{n+1/2} = -\varepsilon^2 \frac{(\kappa_{\mathcal{L}}(\underline{\mathbf{x}}, \underline{\mathbf{w}}))^{n+1} - (\kappa_{\mathcal{L}}(\underline{\mathbf{x}}, \underline{\mathbf{w}}))^n}{\Delta t} + \varepsilon (\eta_{\mathcal{L}}(\underline{\mathbf{x}}, \underline{e}, \underline{\mathbf{w}}))^{n+1/2} \\ + \frac{\varepsilon^2}{\Delta t} \tau_{\mathcal{L}}(\mathbf{x}^n, \mathbf{x}^{n+1}, \mathbf{w}^n, \mathbf{w}^{n+1}), $$

where $\kappa_{\mathcal{L}}$ is given by

$$ \kappa_{\mathcal{L}}(\mathbf{x}, \mathbf{v}) := -\frac{1}{b(\mathbf{x})} \mathcal{L}(\mathbf{x})(\mathbf{v}^{\perp}), $$

whereas the couple $(\eta_{\mathcal{L}}, \tau_{\mathcal{L}})$ is defined as

$$ \eta_{\mathcal{L}}(\mathbf{x}, \underline{e}, \mathbf{v}) := - \left( d_{\mathbf{x}} \left( \frac{\mathcal{L}}{b} \right)(\mathbf{x})(\mathbf{v}) \right) (\mathbf{v}^{\perp}) - \frac{1}{b(\mathbf{x})} \mathcal{L}(\mathbf{x})(\mathbf{E}^{\perp}(\mathbf{x})) \\ + \varepsilon \left( e - \frac{1}{2} \|\mathbf{v}\|^2 \right) \left( d_{\mathbf{x}} \left( \frac{\mathcal{L}}{b} \right)(\mathbf{x}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right)(\mathbf{x}) \right) \right) (\mathbf{v}^{\perp}), $$

and

$$ \tau_{\mathcal{L}}(\mathbf{x}, \mathbf{x}', \mathbf{v}, \mathbf{v}') := - \left( \frac{\mathcal{L}}{b}(\mathbf{x}') - \frac{\mathcal{L}}{b}(\mathbf{x}) - d_{\mathbf{x}} \left( \frac{\mathcal{L}}{b} \right) \left( \frac{\mathbf{x}'+\mathbf{x}}{2} \right) (\mathbf{x}' - \mathbf{x}) \right) \left( \left( \frac{\mathbf{v}'+\mathbf{v}}{2} \right)^{\perp} \right) \\ - \left( \frac{1}{2} \left( \frac{\mathcal{L}}{b}(\mathbf{x}') + \frac{\mathcal{L}}{b}(\mathbf{x}) \right) - \frac{\mathcal{L}}{b} \left( \frac{\mathbf{x}'+\mathbf{x}}{2} \right) \right) ((\mathbf{v}')^{\perp} - \mathbf{v}^{\perp}). $$

**Lemma 4.2.** Consider $\mathcal{B}$ a differentiable function of space with values in bilinear maps from $\mathbb{R}^2 \times \mathbb{R}^2$ to some $\mathbb{R}^n$, $n \in \mathbb{N}^*$. Then solutions to (1.4) satisfy

$$ (\mathcal{B}_{\underline{\mathbf{x}}}(\underline{\mathbf{w}}, \underline{\mathbf{w}}))^{n+1/2} = \left( \frac{1}{2} \|\underline{\mathbf{w}}\|^2 \operatorname{Tr}(\mathcal{B}_{\underline{\mathbf{x}}}) \right)^{n+1/2} - \varepsilon^2 \frac{(\kappa_{\mathcal{B}}(\underline{\mathbf{x}}, \underline{\mathbf{w}}))^{n+1} - (\kappa_{\mathcal{B}}(\underline{\mathbf{x}}, \underline{\mathbf{w}}))^n}{\Delta t} + \varepsilon (\eta_{\mathcal{B}}(\underline{\mathbf{x}}, \underline{e}, \underline{\mathbf{w}}))^{n+1/2} \\ + \frac{\varepsilon^2}{\Delta t} \tau_{\mathcal{B}}(\mathbf{x}^n, \mathbf{x}^{n+1}, \mathbf{w}^n, \mathbf{w}^{n+1}), $$

where

$$ \kappa_{\mathcal{B}}(\mathbf{x}, \mathbf{v}) := -\frac{1}{2} \frac{1}{b(\mathbf{x})} \mathcal{B}(\mathbf{x})(\mathbf{v}, \mathbf{v}^{\perp}), $$

whereas $\eta_{\mathcal{B}}$ is given as

$$ \eta_{\mathcal{B}}(\mathbf{x}, \underline{e}, \mathbf{v}) := -\frac{1}{2} \left( d_{\mathbf{x}} \left( \frac{\mathcal{B}}{b} \right)(\mathbf{x})(\mathbf{v}) \right) (\mathbf{v}, \mathbf{v}^{\perp}) - \frac{1}{2} \frac{1}{b(\mathbf{x})} \left( \mathcal{B}(\mathbf{x})(\mathbf{v}, \mathbf{E}^{\perp}(\mathbf{x})) + \mathcal{B}(\mathbf{x})(\mathbf{E}(\mathbf{x}), \mathbf{v}^{\perp}) \right) \\ + \frac{1}{2} \varepsilon \left( e - \frac{1}{2} \|\mathbf{v}\|^2 \right) \left( d_{\mathbf{x}} \left( \frac{\mathcal{B}}{b} \right)(\mathbf{x}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right)(\mathbf{x}) \right) \right) (\mathbf{v}, \mathbf{v}^{\perp}), $$
and

$$
\begin{align*}
\tau_B(\mathbf{x}, \mathbf{x}', \mathbf{v}, \mathbf{v}') &:= -\frac{1}{2} \left( \frac{\mathcal{B}}{b}(\mathbf{x}') - \frac{\mathcal{B}}{b}(\mathbf{x}) - \mathrm{d}_\mathbf{x} \left( \frac{\mathcal{B}}{b} \right) \left( \frac{\mathbf{x}'+\mathbf{x}}{2} \right) (\mathbf{x}' - \mathbf{x}) \right) \left( \frac{\mathbf{v}'+\mathbf{v}}{2}, \left( \frac{\mathbf{v}'+\mathbf{v}}{2} \right)^\perp \right) \\
&\quad - \frac{1}{2} \left( \frac{1}{2} \left( \frac{\mathcal{B}}{b}(\mathbf{x}') + \frac{\mathcal{B}}{b}(\mathbf{x}) \right) - \frac{\mathcal{B}}{b} \left( \frac{\mathbf{x}'+\mathbf{x}}{2} \right) \right) \left( \frac{\mathbf{v}'+\mathbf{v}}{2}, (\mathbf{v}')^\perp - \mathbf{v}^\perp \right) \\
&\quad - \frac{1}{2} \left( \frac{1}{2} \left( \frac{\mathcal{B}}{b}(\mathbf{x}') + \frac{\mathcal{B}}{b}(\mathbf{x}) \right) - \frac{\mathcal{B}}{b} \left( \frac{\mathbf{x}'+\mathbf{x}}{2} \right) \right) \left( \mathbf{v}' - \mathbf{v}, \left( \frac{\mathbf{v}'+\mathbf{v}}{2} \right)^\perp \right) \\
&\quad - \frac{1}{8} \left( \frac{\mathcal{B}}{b}(\mathbf{x}') - \frac{\mathcal{B}}{b}(\mathbf{x}) \right) (\mathbf{v}' - \mathbf{v}, (\mathbf{v}')^\perp - \mathbf{v}^\perp) .
\end{align*}
$$

Here we apply the foregoing lemmas with exactly the same choices as in the proof of Proposition 2.3, explicitly $\mathcal{L}_1$ constant with value the identity map and $\mathcal{B}_1$ defined by

$$ \mathcal{B}_1(\mathbf{x})(\mathbf{v}_1, \mathbf{v}_2) := - \mathrm{d}_\mathbf{x} \left( \frac{1}{b} \right) (\mathbf{x})(\mathbf{v}_1) \mathbf{v}_2^\perp. $$

This allows us to transform the $\mathbf{x}$ equation of (1.4) into

$$ (4.2) \qquad \frac{(\mathbf{x} + \varepsilon \kappa_1)^{n+1} - (\mathbf{x} + \varepsilon \kappa_1)^n}{\Delta t} = \left( \mathbf{F}_{\underline{\mathbf{x}}} - e \left( \nabla_{\underline{\mathbf{x}}} \left( \frac{1}{b} \right) \right)_{\underline{\mathbf{x}}} \right)^{n+1/2} + \varepsilon (\eta_1^n + \tau_1^n), $$

with for any $n$

$$ \|\kappa_1^n\| \lesssim \|\mathbf{w}^n\| (1 + \varepsilon \|\mathbf{w}^n\|) $$

$$ \|\eta_1^n\| \lesssim \|\underline{\mathbf{w}}^{n+1/2}\| + \|\underline{\mathbf{w}}^{n+1/2}\| (1 + \varepsilon \|\underline{\mathbf{w}}^{n+1/2}\|) (|e^{n+1/2}| + \|\underline{\mathbf{w}}^{n+1/2}\|^2) $$

$$ \|\tau_1^n\| \lesssim \frac{\|\mathbf{x}^{n+1} - \mathbf{x}^n\|}{\Delta t} \left( \|\mathbf{x}^{n+1} - \mathbf{x}^n\|^2 \|\underline{\mathbf{w}}^{n+1/2}\| (1 + \varepsilon \|\underline{\mathbf{w}}^{n+1/2}\|) \right. \\ \left. + \|\mathbf{x}^{n+1} - \mathbf{x}^n\| \|\mathbf{w}^{n+1} - \mathbf{w}^n\| (1 + \varepsilon \|\underline{\mathbf{w}}^{n+1/2}\|) + \varepsilon \|\mathbf{w}^{n+1} - \mathbf{w}^n\|^2 \right). $$

This is sufficient to prove the following discrete counterpart to Proposition 2.3.

**Proposition 4.3.** Let $T > 0$, $M > 0$, $\varepsilon_0 > 0$. There exist $\delta_0 > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, e, \mathbf{w})$ solves (1.4) with initial datum $(\mathbf{x}^0, e^0, \mathbf{w}^0)$ such that $\|\mathbf{v}^0\| \le M$, $|e^0| \le \frac{1}{2}M^2$ and $(\mathbf{y}, g)$ solves (1.5) with initial datum $(\mathbf{y}^0, g^0) = (\mathbf{x}^0, e^0)$ then when $0 \le n \Delta t \le T$,

$$ \|(\mathbf{x}^n, e^n) - (\mathbf{y}^n, g^n)\| \le C \varepsilon. $$

*Proof.* To begin with we invoke Proposition 3.2 to conclude that when $0 \le n \Delta t \le T$ and $\Delta t$ is sufficiently small (independently of $\varepsilon$), $\kappa_1^n$ and $\eta_1^n$ are uniformly bounded whereas

$$ \|\tau_1^n\| \lesssim \min \left( \left\{ \frac{(\Delta t + \varepsilon)^2}{\Delta t} ; \frac{\Delta t}{\varepsilon^2} \right\} \right) \lesssim \min \left( \left\{ \Delta t + \frac{1}{\lambda} ; \lambda \right\} \right) \lesssim 1 + \Delta t, $$

where again $\lambda = \Delta t / \varepsilon^2$.

One may then proceed as in the proof of Proposition 2.3, summing (4.2) instead of integrating (2.4), to deduce that for some constant $K$ (depending on $(T, M, \varepsilon_0)$ and the fields but not on $(\varepsilon, \Delta t)$), provided that $0 \le n \Delta t \le T$ and $\Delta t$ is sufficiently small (independently of $\varepsilon$),

$$ \|\mathbf{x}^n - \mathbf{y}^n\| \le K \Delta t \sum_{\ell=1}^{n} \|\mathbf{x}^{\ell} - \mathbf{y}^{\ell}\| + K \varepsilon, $$
which implies provided that $\Delta t$ is small enough to enforce $K \Delta t \le \frac{1}{2}$,

$$ \|\mathbf{x}^n - \mathbf{y}^n\| \le 2K \Delta t \sum_{\ell=1}^{n-1} \|\mathbf{x}^\ell - \mathbf{y}^\ell\| + 2K \varepsilon, $$

One then concludes with a discrete Grönwall lemma. Explicitly, set

$$ \xi^n := \sum_{\ell=1}^{n-1} \|\mathbf{x}^\ell - \mathbf{y}^\ell\| + \frac{\varepsilon}{\Delta t} $$

so that for any $n$ such that $0 \le n \Delta t \le T$

$$ \xi^{n+1} \le (1 + 2K \Delta t) \xi^n. $$

This implies recursively that for any $n$ such that $0 \le n \Delta t \le T$

$$ \xi^n \le e^{2K n \Delta t} \xi^0 \le e^{2K T} \frac{\varepsilon}{\Delta t}. $$

In turn one deduces that for any $n$ such that $0 \le n \Delta t \le T$

$$ \|\mathbf{x}^n - \mathbf{y}^n\| \le 2K \Delta t \xi^n \le \varepsilon 2K e^{2K T}. $$

Since $e^n - g^n = -(\phi(\mathbf{x}^n) - \phi(\mathbf{y}^n))$, this concludes the proof with $C = 2K e^{2K T} (1 + \|\nabla_{\mathbf{x}}\phi\|_{L^\infty})$. $\square$

## 5. DIRECT CONVERGENCE ANALYSIS

Before diving into the direct convergence analysis a few comments are in order. The general strategy follows the classical path of analyzing separately consistency errors and stability estimates. For background on classical numerical analysis of ordinary differential equations we refer for instance to [14, 15].

On the consistency side the basic observations here are that for any smooth function of time $a$

$$ \frac{a(t + \Delta t) + a(t)}{2} - a\left(t + \frac{1}{2}\Delta t\right) = \frac{1}{8}(\Delta t)^2 \int_{-1}^{1} \left(\frac{d}{dt}\right)^2 a\left(t + \frac{1}{2}(1 + \tau)\Delta t\right) (1 - |\tau|) d\tau $$

$$ \frac{a(t + \Delta t) - a(t)}{\Delta t} - \frac{d a}{d t} \left(t + \frac{1}{2} \Delta t\right) = \frac{1}{16} (\Delta t)^2 \int_{-1}^{1} \left(\frac{d}{dt}\right)^3 a \left(t + \frac{1}{2}(1 + \tau)\Delta t\right) (1 - |\tau|)^2 d\tau. $$

Thus if $a$ is solving

$$ \frac{d a}{d t} = \mathcal{V}(a) $$

for some vector field $\mathcal{V}$ then

$$ \left\| \frac{a(t + \Delta t) - a(t)}{\Delta t} - \mathcal{V} \left( \frac{a(t + \Delta t) + a(t)}{2} \right) \right\| \le \frac{(\Delta t)^2}{8} \left( \left\| \nabla_a \mathcal{V} \right\|_{L^\infty} \left\| \left( \frac{d}{dt} \right)^2 a \right\|_{L^\infty} + \frac{1}{3} \left\| \left( \frac{d}{dt} \right)^3 a \right\|_{L^\infty} \right). $$

This rough classical computation is sufficient to justify that we are indeed studying a second-order scheme and if we were not interested in the dependency on $\varepsilon$ or if we were studying the regime when $\varepsilon$ is bounded away from zero it would be sufficient to conclude.

It is also sufficient to predict what are the best consistency estimates we can expect. Indeed there is indeed no way, even through a near identity change of variable, to improve on the fact that the $k$-th derivative of $\mathbf{v}$ scales like $\varepsilon^{-2k}$ so that the consistency error for $\mathbf{v}$ cannot be better than $(\Delta t)^2 \varepsilon^{-6}$. In turn this implies that the consistency error for a slow variable coupled with $\varepsilon^m \mathbf{v}$ for some $m$ cannot be improved beyond $(\Delta t)^2 \varepsilon^{m-6}$. In particular the consistency error for $\mathbf{x}$ is capped by $(\Delta t)^2 \varepsilon^{-5}$ whereas the one for $\mathbf{x}_{gc}$ cannot beat $(\Delta t)^2 \varepsilon^{-4}$. For a more general discussion we again refer to [10, Appendix A]. In particular, for the sake of comparison note that when [10] analyzes a first-order scheme the corresponding ceilings are respectively $(\Delta t) \varepsilon^{-3}$ and $(\Delta t) \varepsilon^{-2}$.
On the stability side note that the implication is designed to improve how perturbations on fast variables impact both slow and fast variables. The main remaining limitation then arises from how slow variables affect fast variables. The corresponding Lipschitz constant is mostly due to the fact that the amplitude of the magnetic field is non constant and scales as $\epsilon^{-2}$. As a consequence one cannot start the convergence analysis before a level of preparation comparable to the one achieved to prove Proposition 2.5 so as to enforce that the analysis is performed with slow variables and $\epsilon^2 \mathbf{v}$ thus taming the slow-to-fast stiffness. For the sake of comparison we recall that [10] analyzes the case when $b$ is constant so that the slow-to-fast Lipschitz constant scales as $\epsilon^{-1}$ and the analysis may be carried out with slow variables and $\epsilon \mathbf{v}$, requiring only a cheaper preparation, similar to the one to prove Propositions 2.3 and 4.3.

Let us stress that this means that even if one is only interested in $\mathbf{x}$ (and not in $\mathbf{x}_{gc}$) heeding to stability considerations one must nevertheless perform a higher order preparation. This is the reason why we have expounded Proposition 2.5 and its proof even if we omit to state and prove its discrete counterpart.

### 5.1. Stability analysis.
As in Section 2.2 we introduce a guiding center variable

$$ \mathbf{x}_{gc} := \mathbf{x} - \varepsilon \frac{\mathbf{w}^{\perp}}{b(\mathbf{x})} $$

and the associated modified kinetic energy

$$ e_{gc} := e + \frac{\varepsilon}{b(\mathbf{x})} \mathbf{E}^{\perp}(\mathbf{x}) \cdot \mathbf{w}. $$

Note that with $\kappa_e$ as in (2.5)

$$ (5.1) \qquad e_{gc}^{n+1} + \phi(\mathbf{x}_{gc}^{n+1}) + \varepsilon^2 \kappa_e(\mathbf{x}^{n+1}, \mathbf{w}^{n+1}) = e_{gc}^n + \phi(\mathbf{x}_{gc}^n) + \varepsilon^2 \kappa_e(\mathbf{x}^n, \mathbf{w}^n). $$

To proceed we need both a discrete version of Lemma 2.4 and an extension of Lemma 4.1 allowing for linear operators depending on more slow variables. To facilitate the derivation of the latter note that solutions to (1.4) satisfy

$$ \frac{\mathbf{x}_{gc}^{n+1} - \mathbf{x}_{gc}^n}{\Delta t} = - \left( \left( \left( e - \frac{1}{2} \|\mathbf{w}\|^2 \right) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \right)_{\mathbf{x}} \right)^{n+1/2} + (\eta_{\mathcal{L}_1}(\mathbf{x}, e, \mathbf{w}))^{n+1/2} \right. \\ \left. + \frac{\varepsilon}{\Delta t} \tau_{\mathcal{L}_1}(\mathbf{x}^n, \mathbf{x}^{n+1}, \mathbf{w}^n, \mathbf{w}^{n+1}), \right) $$

with notation from Lemma 4.1 and $\mathcal{L}_1 \equiv \text{Id}$.

**Lemma 5.1.** Consider $\mathcal{T}$ a differentiable function of space with values in symmetric trilinear maps from $\mathbb{R}^2 \times \mathbb{R}^2 \times \mathbb{R}^2$ to some $\mathbb{R}^n$, $n \in \mathbb{N}^*$. Then solutions to (1.4) satisfy

$$ (\mathcal{T}_{\underline{\mathbf{x}}}(\underline{\mathbf{w}}, \underline{\mathbf{w}}, \underline{\mathbf{w}}))^{n+1/2} = -\varepsilon^2 \frac{(\kappa \mathcal{T}(\underline{\mathbf{x}}, \underline{\mathbf{w}}))^{n+1} - (\kappa \mathcal{T}(\underline{\mathbf{x}}, \underline{\mathbf{w}}))^n}{\Delta t} + \varepsilon (\eta \mathcal{T}(\underline{\mathbf{x}}, e, \underline{\mathbf{w}}))^{n+1/2} \\ + \frac{\varepsilon^2}{\Delta t} \tau \mathcal{T}(\underline{\mathbf{x}}^n, \underline{\mathbf{x}}^{n+1}, \underline{\mathbf{w}}^n, \underline{\mathbf{w}}^{n+1}), $$

where

$$ \kappa \mathcal{T}(\underline{\mathbf{x}}, \underline{\mathbf{v}}) := - \frac{1}{b(\underline{\mathbf{x}})} \left( \mathcal{T}(\underline{\mathbf{x}})(\underline{\mathbf{v}}, \underline{\mathbf{v}}, \underline{\mathbf{v}}^{\perp}) + \frac{2}{3} \mathcal{T}(\underline{\mathbf{x}})(\underline{\mathbf{v}}^{\perp}, \underline{\mathbf{v}}^{\perp}, \underline{\mathbf{v}}^{\perp}) \right), $$
and

$$
\begin{align*}
\eta \mathcal{T}(\mathbf{x}, e, \mathbf{v}) &:= - \left( d_{\mathbf{x}} \left( \frac{\mathcal{T}}{b} \right)(\mathbf{x})(\mathbf{v}) \right) (\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) - \frac{2}{3} \left( d_{\mathbf{x}} \left( \frac{\mathcal{T}}{b} \right)(\mathbf{x})(\mathbf{v}) \right) (\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp}) \\
&\quad - \frac{1}{b(\mathbf{x})} \left( \mathcal{T}(\mathbf{x})(\mathbf{v}, \mathbf{v}, \mathbf{E}^{\perp}(\mathbf{x})) + 2\mathcal{T}(\mathbf{x})(\mathbf{E}(\mathbf{x}), \mathbf{v}, \mathbf{v}^{\perp}) + 2\mathcal{T}(\mathbf{x})(\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{E}^{\perp}(\mathbf{x})) \right) \\
&\quad + \varepsilon \left( e - \frac{1}{2} \|\mathbf{v}\|^2 \right) \left( d_{\mathbf{x}} \left( \frac{\mathcal{T}}{b} \right)(\mathbf{x}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right)(\mathbf{x}) \right) \right) (\mathbf{v}, \mathbf{v}, \mathbf{v}^{\perp}) \\
&\quad + \frac{2}{3} \varepsilon \left( e - \frac{1}{2} \|\mathbf{v}\|^2 \right) \left( d_{\mathbf{x}} \left( \frac{\mathcal{T}}{b} \right)(\mathbf{x}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right)(\mathbf{x}) \right) \right) (\mathbf{v}^{\perp}, \mathbf{v}^{\perp}, \mathbf{v}^{\perp}),
\end{align*}
$$

$$
\begin{align*}
\tau \mathcal{T}(\mathbf{x}, \mathbf{x}', \mathbf{v}, \mathbf{v}') &:= - \left( \frac{\mathcal{T}}{b}(\mathbf{x}') - \frac{\mathcal{T}}{b}(\mathbf{x}) - d_{\mathbf{x}} \left( \frac{\mathcal{T}}{b} \right) \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) (\mathbf{x}' - \mathbf{x}) \right) \left( \frac{\mathbf{v}' + \mathbf{v}}{2}, \frac{\mathbf{v}' + \mathbf{v}}{2}, \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp} \right) \\
&\quad - \frac{2}{3} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') - \frac{\mathcal{T}}{b}(\mathbf{x}) - d_{\mathbf{x}} \left( \frac{\mathcal{T}}{b} \right) \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) (\mathbf{x}' - \mathbf{x}) \right) \left( \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp}, \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp}, \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp} \right) \\
&\quad - \left( \frac{1}{2} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') + \frac{\mathcal{T}}{b}(\mathbf{x}) \right) - \frac{\mathcal{T}}{b} \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \right) \left( \frac{\mathbf{v}' + \mathbf{v}}{2}, \frac{\mathbf{v}' + \mathbf{v}}{2}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp} \right) \\
&\quad - 2 \left( \frac{1}{2} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') + \frac{\mathcal{T}}{b}(\mathbf{x}) \right) - \frac{\mathcal{T}}{b} \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \right) \left( \mathbf{v}' - \mathbf{v}, \frac{\mathbf{v}' + \mathbf{v}}{2}, \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp} \right) \\
&\quad - 2 \left( \frac{1}{2} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') + \frac{\mathcal{T}}{b}(\mathbf{x}) \right) - \frac{\mathcal{T}}{b} \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \right) \left( \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp}, \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp} \right) \\
&\quad - \frac{1}{8} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') + \frac{\mathcal{T}}{b}(\mathbf{x}) \right) \left( \mathbf{v}' - \mathbf{v}, \mathbf{v}' - \mathbf{v}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp} \right) \\
&\quad - \frac{1}{12} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') + \frac{\mathcal{T}}{b}(\mathbf{x}) \right) \left( (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp} \right) \\
&\quad - \frac{1}{4} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') - \frac{\mathcal{T}}{b}(\mathbf{x}) \right) \left( \mathbf{v}' - \mathbf{v}, \mathbf{v}' - \mathbf{v}, \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp} \right) \\
&\quad - \frac{1}{2} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') - \frac{\mathcal{T}}{b}(\mathbf{x}) \right) \left( \frac{\mathbf{v}' + \mathbf{v}}{2}, \mathbf{v}' - \mathbf{v}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp} \right) \\
&\quad - \frac{1}{2} \left( \frac{\mathcal{T}}{b}(\mathbf{x}') - \frac{\mathcal{T}}{b}(\mathbf{x}) \right) \left( \left( \frac{\mathbf{v}' + \mathbf{v}}{2} \right)^{\perp}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp}, (\mathbf{v}')^{\perp} - \mathbf{v}^{\perp} \right).
\end{align*}
$$

Again the proof of Lemma 5.1 is a straightforward — but lengthy ! — computation, that we mostly omit. To provide some elements of it, let us only stress that if $\mathbf{T}_0$ is a fixed symmetric trilinear map then

$$
\begin{align*}
(\mathbf{T}_0(a, a, b))^{n+1} - (\mathbf{T}_0(a, a, b))^n
  &= \mathbf{T}_0(\underline{a}^{n+1/2}, \underline{a}^{n+1/2}, b^{n+1} - b^n) + \frac{1}{4} \mathbf{T}_0(a^{n+1} - a^n, a^{n+1} - a^n, b^{n+1} - b^n) \\
  &\quad + 2 \mathbf{T}_0(a^{n+1} - a^n, \underline{a}^{n+1/2}, b^{n+1/2})
\end{align*}
$$
and

$$ \underline{\mathbf{T}_0(a, a, b)}^{n+1/2} = \mathbf{T}_0(\underline{a}^{n+1/2}, \underline{a}^{n+1/2}, \underline{b}^{n+1/2}) + \frac{1}{4} \mathbf{T}_0(a^{n+1} - a^n, a^{n+1} - a^n, \underline{b}^{n+1/2}) \\ + \frac{1}{2} \mathbf{T}_0(\underline{a}^{n+1/2}, a^{n+1} - a^n, \underline{b}^{n+1} - b^n) $$

thus in particular

$$ (\mathbf{T}_0(a, a, a))^{n+1} - (\mathbf{T}_0(a, a, a))^n \\ = 3 \mathbf{T}_0(\underline{a}^{n+1/2}, \underline{a}^{n+1/2}, a^{n+1} - a^n) + \frac{1}{4} \mathbf{T}_0(a^{n+1} - a^n, a^{n+1} - a^n, a^{n+1} - a^n) $$

and

$$ \underline{\mathbf{T}_0(a, a, a)}^{n+1/2} = \mathbf{T}_0(\underline{a}^{n+1/2}, \underline{a}^{n+1/2}, \underline{a}^{n+1/2}) + \frac{3}{4} \mathbf{T}_0(\underline{a}^{n+1/2}, a^{n+1} - a^n, \underline{a}^{n+1} - a^n). $$

**Lemma 5.2.** Consider $\mathcal{L}$ a differentiable function of space-guiding center-kinetic energy with values in linear maps from $\mathbb{R}^2$ to some $\mathbb{R}^n$, $n \in \mathbb{N}^*$. Then solutions to (1.4) satisfy

$$ (\mathcal{L}_{\underline{x},\underline{x}_{gc},\underline{e}}(\underline{\mathbf{w}}))^{n+1/2} = -\varepsilon^2 \frac{(\kappa\mathcal{L}(\underline{x}, \underline{x}_{gc}, e, \underline{\mathbf{w}}))^{n+1} - (\kappa\mathcal{L}(\underline{x}, \underline{x}_{gc}, e, \underline{\mathbf{w}}))^n}{\Delta t} + \varepsilon (\eta\mathcal{L}(\underline{x}, \underline{x}_{gc}, \underline{e}, \underline{\mathbf{w}}))^{n+1/2} \\ + \frac{\varepsilon^2}{\Delta t} \tau_{\mathcal{L}}(\underline{x}^n, \underline{x}^{n+1}, \underline{x}_{gc}^n, \underline{x}_{gc}^{n+1}, e^n, e^{n+1}, \underline{\mathbf{w}}^n, \underline{\mathbf{w}}^{n+1}), $$

where

$$ \kappa_{\mathcal{L}}(\underline{x}, \underline{y}, e, \underline{v}) := -\frac{1}{b(\underline{x})} \mathcal{L}(\underline{x}, \underline{y}, e)(\underline{v}^{\perp}), $$

and

$$ \eta_{\mathcal{L}}(\underline{x}, \underline{y}, e, \underline{v}) \\ := - \left( d_{\underline{x}} \left( \frac{\mathcal{L}}{b} \right) (\underline{x}, \underline{y}, e)(\underline{v}) \right) (\underline{v}^{\perp}) - \frac{1}{b(\underline{x})} \mathcal{L}(\underline{x}, \underline{y}, e)(\underline{E}^{\perp}(\underline{x})) \\ + \varepsilon \left( e - \frac{1}{2} \|\underline{v}\|^2 \right) \left( d_{\underline{x}} \left( \frac{\mathcal{L}}{b} \right) (\underline{x}, \underline{y}, e) \left( \nabla_{\underline{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{x}) \right) \right) (\underline{v}^{\perp}) \\ - \frac{1}{b(\underline{x})} \left( d_e \mathcal{L}(\underline{x}, \underline{y}, e)(\underline{E}(\underline{x}) \cdot \underline{v}) \right) (\underline{v}^{\perp}) \\ + \frac{\varepsilon}{b(\underline{x})} \left( e - \frac{1}{2} \|\underline{v}\|^2 \right) \left( d_e \mathcal{L}(\underline{x}, \underline{y}, e) \left( \underline{E}(\underline{x}) \cdot \nabla_{\underline{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{x}) \right) \right) (\underline{v}^{\perp}) \\ - \varepsilon \left( d_{\underline{y}} \left( \frac{\mathcal{L}}{b} \right) (\underline{x}, \underline{y}, e) \left( - \left( e - \frac{1}{2} \|\underline{v}\|^2 \right) \nabla_{\underline{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{x}) + \eta_{\mathcal{L}_1}(\underline{x}, e, \underline{v}) \right) \right) (\underline{v}^{\perp}), $$
$$
\begin{align*}
& \tau_{\mathcal{L}}(\mathbf{x}, \mathbf{x}', \mathbf{y}, \mathbf{y}', e, e', \mathbf{v}, \mathbf{v}') \\
&:= - \left( \frac{\mathcal{L}}{b}(\mathbf{x}', \mathbf{y}', e') - \frac{\mathcal{L}}{b}(\mathbf{x}, \mathbf{y}, e) \right. \\
& \qquad \left. - d_{(\mathbf{x}, \mathbf{y}, e)} \left( \frac{\mathcal{L}}{b} \right) \left( \frac{\mathbf{x}'+\mathbf{x}}{2}, \frac{\mathbf{y}'+\mathbf{y}}{2}, \frac{e'+e}{2} \right) ((\mathbf{x}' - \mathbf{x}, \mathbf{y}' - \mathbf{y}, e' - e)) \right) \left( \left( \frac{\mathbf{v}'+\mathbf{v}}{2} \right)^{\perp} \right) \\
& \quad - \left( \frac{1}{2} \left( \frac{\mathcal{L}}{b}(\mathbf{x}', \mathbf{y}', e') + \frac{\mathcal{L}}{b}(\mathbf{x}, \mathbf{y}, e) \right) - \frac{\mathcal{L}}{b} \left( \frac{\mathbf{x}'+\mathbf{x}}{2}, \frac{\mathbf{y}'+\mathbf{y}}{2}, \frac{e'+e}{2} \right) \right) ((\mathbf{v}')^{\perp} - \mathbf{v}^{\perp}) \\
& \quad + \left( d_e \left( \frac{\mathcal{L}}{b} \right) \left( \frac{\mathbf{x}'+\mathbf{x}}{2}, \frac{\mathbf{y}'+\mathbf{y}}{2}, \frac{e'+e}{2} \right) \left( \phi(\mathbf{x}') - \phi(\mathbf{x}) - d_{\mathbf{x}} \phi \left( \frac{\mathbf{x}'+\mathbf{x}}{2} \right) (\mathbf{x}' - \mathbf{x}) \right) \right) \left( \left( \frac{\mathbf{v}'+\mathbf{v}}{2} \right)^{\perp} \right) \\
& \quad - \left( d_{\mathbf{y}} \left( \frac{\mathcal{L}}{b} \right) \left( \frac{\mathbf{x}'+\mathbf{x}}{2}, \frac{\mathbf{y}'+\mathbf{y}}{2}, \frac{e'+e}{2} \right) (\tau_{\mathcal{L}_1}(\mathbf{x}, \mathbf{x}', \mathbf{v}, \mathbf{v}')) \right) \left( \left( \frac{\mathbf{v}'+\mathbf{v}}{2} \right)^{\perp} \right).
\end{align*}
$$

The proof of Lemma 5.2 is a straightforward variation on the proof of Lemma 4.1.

Then we apply Lemma 5.1 to $\mathcal{T}_1$ and Lemma 5.2 to $\mathcal{L}_2$ with $\mathcal{T}_1$ and $\mathcal{L}_2$ defined as in Subsection 2.2 so as to modify (4.2). With the consistency analysis in mind, we want to express the outcome, at leading order, as a scheme for (2.6). To do so we need to compare (4.2) with (2.4) and keep track of differences between Lemmas 5.1 and 5.2 and Lemmas 2.4 and 2.1.

To begin with note that with $\kappa_1$ and $\eta_1$ as in (2.4), equation (4.2) takes the form

$$
\frac{(\mathbf{x} + \varepsilon \kappa_1(\mathbf{x}, \mathbf{w}))^{n+1} - (\mathbf{x} + \varepsilon \kappa_1(\mathbf{x}, \mathbf{w}))^n}{\Delta t} = \left( \mathbf{F}_{\underline{\mathbf{x}}} - e \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \right)_{\underline{\mathbf{x}}} + \varepsilon \eta_1(\underline{\mathbf{x}}, \mathbf{w}) \right)^{n+1/2} + \varepsilon (\tau_1^n + \tilde{\tau}_1^n),
$$

with $\tau_1^n$ as in (4.2) and

$$
\tilde{\tau}_1^n := \frac{1}{2} \left( \left( e - \frac{1}{2} \|\underline{\mathbf{w}}\|^2 \right) \left( d_{\underline{\mathbf{x}}} \left( \frac{\mathcal{B}_1}{b} \right) (\underline{\mathbf{x}}) \left( \nabla_{\underline{\mathbf{x}}}^{\perp} \left( \frac{1}{b} \right) (\underline{\mathbf{x}}) \right) \right) (\underline{\mathbf{w}}, \underline{\mathbf{w}}^{\perp}) \right)^{n+1/2}.
$$

Then with $\kappa_2$ and $\eta_3$ as in (2.6), equation (4.2) is transformed into

$$
(5.2) \quad \frac{(\mathbf{x}_{gc} + \varepsilon^2 \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{w}))^{n+1} - (\mathbf{x}_{gc} + \varepsilon^2 \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{w}))^n}{\Delta t} \\
= \left( \mathbf{F}_{\underline{\mathbf{x}_{gc}}} - e_{gc} \left( \nabla_{\underline{\mathbf{x}}}^{\perp} \left( \frac{1}{b} \right) \right)_{\underline{\mathbf{x}_{gc}}} + \varepsilon^2 \eta_2(\underline{\mathbf{x}}, \mathbf{x}_{gc}, e, \underline{\mathbf{w}}) \right)^{n+1/2} + \varepsilon (\tau_1^n + \tilde{\tau}_1^n) + \varepsilon^2 (\tau_2^n + \tilde{\tau}_2^n),
$$

with

$$
\tau_2^n := \frac{\varepsilon}{\Delta t} \tau_{\mathcal{L}_2}(\mathbf{x}^n, \mathbf{x}^{n+1}, \mathbf{x}_{gc}^n, \mathbf{x}_{gc}^{n+1}, e^n, e^{n+1}, \mathbf{w}^n, \mathbf{w}^{n+1}) + \frac{\varepsilon}{\Delta t} \tau_{\mathcal{T}_1}(\mathbf{x}^n, \mathbf{x}^{n+1}, \mathbf{w}^n, \mathbf{w}^{n+1})
$$
and

$$
\begin{align*}
\tilde{\tau}_2^n := & \left( \left( \underline{\mathbf{e}} - \frac{1}{2} \|\mathbf{w}\|^2 \right) \left( d_x \left( \frac{\mathcal{T}_1}{b} \right) (\underline{\mathbf{x}}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{\mathbf{x}}) \right) \right) (\underline{\mathbf{w}}, \underline{\mathbf{w}}, \underline{\mathbf{w}}^{\perp}) \right)^{n+1/2} \\
& + \frac{2}{3} \left( \left( \underline{\mathbf{e}} - \frac{1}{2} \|\mathbf{w}\|^2 \right) \left( d_x \left( \frac{\mathcal{T}_1}{b} \right) (\mathbf{x}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\mathbf{x}) \right) \right) (\underline{\mathbf{w}}^{\perp}, \underline{\mathbf{w}}^{\perp}, \underline{\mathbf{w}}^{\perp}) \right)^{n+1/2} \\
& + \left( \left( \underline{\mathbf{e}} - \frac{1}{2} \|\mathbf{w}\|^2 \right) \left( \left( d_x \left( \frac{\mathcal{L}_2}{b} \right) + d_y \left( \frac{\mathcal{L}_2}{b} \right) \right) (\underline{\mathbf{x}}, \underline{\mathbf{x}}_{gc}, \underline{e}) \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{\mathbf{x}}) \right) (\underline{\mathbf{w}}^{\perp}) \right) \right)^{n+1/2} \\
& + \left( \left( \underline{\mathbf{e}} - \frac{1}{2} \|\mathbf{w}\|^2 \right) \frac{1}{b(\underline{\mathbf{x}})} \left( d_e \mathcal{L}_2(\underline{\mathbf{x}}, \underline{\mathbf{x}}_{gc}, \underline{e}) \left( \mathbf{E}(\underline{\mathbf{x}}) \cdot \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) (\underline{\mathbf{x}}) \right) (\underline{\mathbf{w}}^{\perp}) \right) \right)^{n+1/2} .
\end{align*}
$$

Consistently, to any sequence $(\mathbf{x}, e, \mathbf{w})$, (not necessarily solving (1.4)) we associate $\mathbf{x}_{gc}$, $e_{gc}$ through

$$
\mathbf{x}_{gc} := \mathbf{x} - \varepsilon \frac{\mathbf{w}^{\perp}}{b(\mathbf{x})}, \qquad e_{gc} := e - \frac{\varepsilon}{b(\mathbf{x})} \mathbf{E}^{\perp}(\mathbf{x}) \cdot \mathbf{w},
$$

and set

$$
\begin{align*}
\mathrm{Res}_x^n(\mathbf{x}, e, \mathbf{w}) := & \frac{(\mathbf{x}_{gc} + \varepsilon^2 \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{w}))^{n+1} - (\mathbf{x}_{gc} + \varepsilon \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{w}))^n}{\Delta t} \\
& - \left( \mathbf{F}_{\underline{\mathbf{x}}_{gc}} - \underline{e}_{gc} \left( \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \right)_{\underline{\mathbf{x}}_{gc}} + \varepsilon^2 \eta_2(\underline{\mathbf{x}}, \mathbf{x}_{gc}, e, \underline{\mathbf{w}}) \right)^{n+1/2}
\end{align*}
$$

$$
\mathrm{Res}_e^n(\mathbf{x}, e, \mathbf{w}) := \frac{(e_{gc} + \phi(\mathbf{x}_{gc}) + \varepsilon^2 \kappa_e(\mathbf{x}, \mathbf{w}))^{n+1} - (e_{gc} + \phi(\mathbf{x}_{gc}) + \varepsilon^2 \kappa_e(\mathbf{x}, \mathbf{w}))^n}{\Delta t}
$$

and

$$
\mathrm{Res}_w^n(\mathbf{x}, e, \mathbf{w}) := \frac{\mathbf{w}^{n+1} - \mathbf{w}^n}{\Delta t} - \left( \varepsilon^{-1}\mathbf{E}_{\underline{\mathbf{x}}} - \varepsilon^{-2}b_{\underline{\mathbf{x}}}\underline{\mathbf{w}}^{\perp} \right)^{n+1/2}.
$$

**Proposition 5.3.** Let $T > 0$, $M' > 0$. There exist $\varepsilon_0 > 0$, $\delta_0 > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, e, \mathbf{w})$ and $(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})$ are such that when $0 \le n \Delta t \le T$,

$$
\|\mathbf{w}^n\| \le M', \quad |e^n| \le \frac{1}{2}(M')^2, \quad \|\tilde{\mathbf{w}}^n\| \le M', \quad |\tilde{e}^n| \le \frac{1}{2}(M')^2,
$$

then when $0 \le n \Delta t \le T$,

$$
\begin{align*}
& \|(\mathbf{x}_{gc}^n, e_{gc}^n) - (\tilde{\mathbf{x}}_{gc}^n, \tilde{e}_{gc}^n)\| + \varepsilon \|(\mathbf{x}^n, e^n) - (\tilde{\mathbf{x}}^n, \tilde{e}^n)\| + \varepsilon^2 \|\mathbf{w}^n - \tilde{\mathbf{w}}^n\| \\
& \le C \left( \|(\mathbf{x}_{gc}^0, e_{gc}^0) - (\tilde{\mathbf{x}}_{gc}^0, \tilde{e}_{gc}^0)\| + \varepsilon^2 \|\mathbf{w}^0 - \tilde{\mathbf{w}}^0\| \right) \\
& \quad + C \max_{0 \le k \le n} \left( \|\mathrm{Res}_x^k(\mathbf{x}, e, \mathbf{w})\| + \|\mathrm{Res}_x^k(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})\| + |\mathrm{Res}_e^k(\mathbf{x}, e, \mathbf{w})| + |\mathrm{Res}_e^k(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})| \right) \\
& \quad + C \varepsilon^2 \max_{0 \le k \le n} \left( \|\mathrm{Res}_w^k(\mathbf{x}, e, \mathbf{w})\| + \varepsilon^2 \|\mathrm{Res}_w^k(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})\| \right).
\end{align*}
$$

*Proof.* Note that when **w** is restricted to a bounded set the map $(\mathbf{x}, e, \mathbf{w}) \mapsto (\mathbf{x}_{gc} - \mathbf{x}, e_{gc} - e)$ has Lipschitz bound $\mathcal{O}(\varepsilon)$ so that when $\varepsilon$ is small enough (depending on the bound on $\mathbf{w}$) by the Banach fixed point theorem one may invert the relation and see $\mathbf{x}, e$ as functions of $(\mathbf{x}_{gc}, e_{gc}, \mathbf{w})$ with Lipschitz dependence encoded by

$$
\|(\mathbf{x}, e) - (\tilde{\mathbf{x}}, \tilde{e})\| \lesssim \|(\mathbf{x}_{gc}, e_{gc}) - (\tilde{\mathbf{x}}_{gc}, \tilde{e}_{gc})\| + \varepsilon \|\mathbf{w} - \tilde{\mathbf{w}}\|.
$$

In particular it is sufficient to estimate $\|(\mathbf{x}_{gc}^n, e_{gc}^n) - (\tilde{\mathbf{x}}_{gc}^n, \tilde{e}_{gc}^n)\| + \varepsilon^2 \|\mathbf{w}^n - \tilde{\mathbf{w}}^n\|$.
Likewise, seeing $\mathbf{x}, e$ as functions of $(\mathbf{x}_{gc}, e_{gc}, \mathbf{w})$ and restricting to bounded sets of $(e_{gc}, \mathbf{w})$ the map $(\mathbf{x}_{gc}, e_{gc}, \mathbf{w}) \mapsto \varepsilon^2(\kappa_2(\mathbf{x}_{gc}, e, \mathbf{w}), \kappa_e(\mathbf{x}, \mathbf{w}))$ has Lipschitz bound $\mathcal{O}(\varepsilon^2)$ so that the relations

$$ \mathbf{y}_{gc} := \mathbf{x}_{gc} + \varepsilon^2 \kappa_2(\mathbf{x}, \mathbf{x}_{gc}, e, \mathbf{w}), $$

$$ g_{gc} := e_{gc} + \varepsilon^2 \kappa_e(\mathbf{x}, e), $$

may be inverted with Lipschitz dependence

$$ \|(\mathbf{x}_{gc}, e_{gc}) - (\tilde{\mathbf{x}}_{gc}, \tilde{e}_{gc})\| \lesssim \|(\mathbf{y}_{gc}, g_{gc}) - (\tilde{\mathbf{y}}_{gc}, \tilde{g}_{gc})\| + \varepsilon^2 \|\mathbf{w} - \tilde{\mathbf{w}}\|. $$

In particular it is sufficient to estimate $\|(\mathbf{y}_{gc}^n, g_{gc}^n) - (\tilde{\mathbf{y}}_{gc}^n, \tilde{g}_{gc}^n)\| + \varepsilon^2 \|\mathbf{w}^n - \tilde{\mathbf{w}}^n\|$.

With this in mind one may interpret Res terms as forcing consistency errors in schemes for $(\mathbf{y}_{gc}, g_{gc}, \mathbf{w})$ and $(\tilde{\mathbf{y}}_{gc}, \tilde{g}_{gc}, \tilde{\mathbf{w}})$ so that heeding to the $\varepsilon$-dependence of equations for some constant $K$ (independent of $(\varepsilon, \Delta t)$), when $0 \le n \Delta t \le T$

$$ \begin{aligned} & \|(\mathbf{y}_{gc}^{n+1}, g_{gc}^{n+1}) - (\tilde{\mathbf{y}}_{gc}^{n+1}, \tilde{g}_{gc}^{n+1})\| + \varepsilon^2 \|\mathbf{w}^{n+1} - \tilde{\mathbf{w}}^{n+1}\| \\ & \le (1 + K\Delta t) \left( \| (\mathbf{y}_{gc}^n, g_{gc}^n) - (\tilde{\mathbf{y}}_{gc}^n, \tilde{g}_{gc}^n)\| + \varepsilon^2 \|\mathbf{w}^n - \tilde{\mathbf{w}}^n\| \right) \\ & \quad + K\Delta t \left( \| (\mathbf{y}_{gc}^{n+1}, g_{gc}^{n+1}) - (\tilde{\mathbf{y}}_{gc}^{n+1}, \tilde{g}_{gc}^{n+1})\| + \varepsilon^2 \|\mathbf{w}^{n+1} - \tilde{\mathbf{w}}^{n+1}\| \right) \\ & \quad + K\Delta t \left( \| \mathrm{Res}_x^n(\mathbf{x}, e, \mathbf{w})\| + \| \mathrm{Res}_x^n(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})\| + |\mathrm{Res}_e^n(\mathbf{x}, e, \mathbf{w})| + |\mathrm{Res}_e^n(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})| \right) \\ & \quad + K\Delta t \varepsilon^2 \left( \| \mathrm{Res}_w^n(\mathbf{x}, e, \mathbf{w})\| + \varepsilon^2 \| \mathrm{Res}_w^n(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})\| \right). \end{aligned} $$

One then concludes as in the proof of Proposition 4.3 by taking $\Delta t$ sufficiently small and invoking a discrete Grönwall lemma. □

**5.2. Consistency analysis.** The goal is obviously to apply Proposition 5.3 to compare a solution to (1.4) with the time discretization of a solution to (1.1). On the latter side the discussion at the beginning of the current section already proves the following lemma.

**Lemma 5.4.** Let $T > 0$, $M > 0$, $\varepsilon_0 > 0$ and $\delta_0 > 0$. There exists $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, \mathbf{v})$ solves (1.1) with initial datum $(\mathbf{x}^0, \mathbf{v}^0)$ such that $\|\mathbf{v}^0\| \le M$ and $(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})$ is defined by

$$ \tilde{\mathbf{x}}^n := \mathbf{x}(n \Delta t), \qquad \tilde{e}^n := \frac{1}{2} \|\mathbf{v}(n \Delta t)\|^2, \qquad \tilde{\mathbf{w}}^n := \mathbf{v}(n \Delta t), $$

then for any $n$ such that $0 \le n \Delta t \le T$,

$$ \| \mathrm{Res}_x^n(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}}) \| + |\mathrm{Res}_e^n(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})| + \varepsilon^2 \| \mathrm{Res}_w^n(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}}) \| \le C \frac{(\Delta t)^2}{\varepsilon^4}. $$

The following lemma completes the estimates of Proposition 3.2 hence providing estimates necessary to apply Proposition 5.3.

**Lemma 5.5.** Let $T > 0$, $M > 0$ and $\varepsilon_0 > 0$. There exist $\delta_0 > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, e, \mathbf{w})$ solves (1.4) with initial datum $(\mathbf{x}^0, e^0, \mathbf{w}^0)$ such that $\|\mathbf{w}^0\| \le M$, $e^0 = \frac{1}{2} \|\mathbf{w}^0\|^2$ then when $0 \le n \Delta t \le T$,

$$ \left| e^{n+1/2} - \frac{1}{2} \|\underline{\mathbf{w}}^{n+1/2}\|^2 \right| \le C \frac{(\Delta t)^2}{\varepsilon^4}. $$

*Proof.* To begin with, note that

$$ e^{n+1/2} - \frac{1}{2} \|\underline{\mathbf{w}}^{n+1/2}\|^2 = \left( e - \frac{1}{2} \|\underline{\mathbf{w}}\|^2 \right)^{n+1/2} + \frac{1}{8} \|\mathbf{w}^{n+1} - \mathbf{w}^n\|^2 $$

so that in view of $\|\mathbf{w}^{n+1} - \mathbf{w}^n\|$ estimates from Proposition 3.2 it is sufficient to bound $e - \frac{1}{2} \|\mathbf{w}\|^2$. Now recall that by taking the scalar product of the $\mathbf{w}$ equation of (1.4) with $\underline{\mathbf{w}}^{n+1/2}$ one obtains

$$ \|\mathbf{w}^{n+1}\|^2 - \|\mathbf{w}^n\|^2 = 2 \Delta t \mathbf{E}_{\underline{\mathbf{x}}}^{n+1/2} \cdot \frac{\mathbf{w}^{n+1/2}}{\varepsilon}, $$
so that

$$ \left(e - \frac{1}{2}\|\mathbf{w}\|^2\right)^{n+1} - \left(e - \frac{1}{2}\|\mathbf{w}\|^2\right)^n = -\Delta t \left(\phi(\mathbf{x}^{n+1}) - \phi(\mathbf{x}^n) - d_{\mathbf{x}} \phi(\mathbf{x}^{n+1/2})(\mathbf{x}^{n+1} - \mathbf{x}^n)\right) \\ - \Delta t \left(\left(e - \frac{1}{2}\|\mathbf{w}\|^2\right) \mathbf{E}_{\underline{\mathbf{x}}} \cdot \nabla_{\underline{\mathbf{x}}}^{\perp} \left(\frac{1}{b}\right)(\underline{\mathbf{x}})\right)^{n+1/2} $$

hence for some constant $K$ (independent of $(\varepsilon, \Delta t)$)

$$ \left| \left(e - \frac{1}{2}\|\mathbf{w}\|^2\right)^{n+1} - \left(e - \frac{1}{2}\|\mathbf{w}\|^2\right)^n \right| \le K \Delta t \left| \left(e - \frac{1}{2}\|\mathbf{w}\|^2\right)^{n+1/2} \right| + K \Delta t \left( \|\mathbf{x}^{n+1} - \mathbf{x}^n\|^2 + \|\mathbf{w}^{n+1} - \mathbf{w}^n\|^2 \right). $$

With estimates from Proposition 3.2 one then concludes as in the proof of Proposition 4.3 by taking $\Delta t$ sufficiently small and invoking a discrete Grönwall lemma. $\square$

**Corollary 5.6.** Let $T > 0$, $M > 0$, $\varepsilon_0 > 0$. There exist $\delta_0 > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, e, \mathbf{w})$ solves (1.4) with initial datum $(\mathbf{x}^0, e^0, \mathbf{w}^0)$ such that $\|\mathbf{w}^0\| \le M$, $e^0 = \frac{1}{2}\|\mathbf{w}^0\|^2$ then when $0 \le n \Delta t \le T$,

$$ \| \mathrm{Res}_x^n(\mathbf{x}, e, \mathbf{w}) \| + |\mathrm{Res}_e^n(\mathbf{x}, e, \mathbf{w})| + \varepsilon^2 \| \mathrm{Res}_w^n(\mathbf{x}, e, \mathbf{w}) \| \le C \frac{(\Delta t)^2}{\varepsilon^4}. $$

### 5.3. Convergence.

**Theorem 5.7.** Let $T > 0$, $M > 0$ and $\varepsilon_0 > 0$. There exist $\delta_0 > 0$ and $C > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, \mathbf{v})$ solves (1.1) with initial datum $(\mathbf{x}^0, \mathbf{v}^0)$ such that $\|\mathbf{v}^0\| \le M$ and $(\tilde{\mathbf{x}}, \tilde{e}, \tilde{\mathbf{w}})$ is defined by

$$ \tilde{\mathbf{x}}^n := \mathbf{x}(n \Delta t), \qquad \tilde{e}^n := \frac{1}{2} \|\mathbf{v}(n \Delta t)\|^2, \qquad \tilde{\mathbf{w}}^n := \mathbf{v}(n \Delta t), $$

and $(\mathbf{x}, e, \mathbf{w})$ solves (1.4) with initial datum $(\mathbf{x}^0, e^0, \mathbf{w}^0)$ such that $\mathbf{w}^0 = \mathbf{v}^0$, $e^0 = \frac{1}{2}\|\mathbf{v}^0\|^2$, then for any $n$ such that $0 \le n \Delta t \le T$,

$$ \|(\mathbf{x}_{\mathrm{gc}}^n, e_{\mathrm{gc}}^n) - (\tilde{\mathbf{x}}_{\mathrm{gc}}^n, \tilde{e}_{\mathrm{gc}}^n)\| + \varepsilon \|(\mathbf{x}^n, e^n) - (\tilde{\mathbf{x}}^n, \tilde{e}^n)\| + \varepsilon^2 \|\mathbf{w}^n - \tilde{\mathbf{w}}^n\| \le C \frac{(\Delta t)^2}{\varepsilon^4}. $$

*Proof.* Note that Section 2 and Proposition 3.2 provide the bounds necessary to apply Proposition 5.3. Therefore by combining Proposition 5.3 with Lemma 5.4 and Corollary 5.6 provides the theorem under the extra assumption (needed to apply Proposition 5.3) that $0 < \varepsilon \le \varepsilon'_0$ for some sufficiently small $\varepsilon'_0$.

To complete the proof it is sufficient to observe that when $\varepsilon'_0 > 0$ is held fixed and $\varepsilon'_0 \le \varepsilon \le \varepsilon_0$ standard direct convergence analysis in original variables $(\mathbf{x}, e, \mathbf{w})$ provides the missing bounds (with constants blowing up dramatically in the limit $\varepsilon'_0 \to 0$). $\square$

## 6. PROOF OF THE MAIN THEOREM, THEOREM 1.1

In this short section we explain how to combine the previous elements so as to conclude the proof of Theorem 1.1. Half of the conclusion is the bound

$$ \|(\mathbf{x}^n, e^n) - (\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2) (n \Delta t)\| + \varepsilon \|\mathbf{w}^n - \mathbf{v}(n \Delta t)\| \le C \frac{(\Delta t)^2}{\varepsilon^5}, $$

which is a direct consequence of Theorem 5.7.

The second half is

$$ \|(\mathbf{x}^n, e^n) - (\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2) (n \Delta t)\| + \varepsilon \|\mathbf{w}^n - \mathbf{v}(n \Delta t)\| \le C (\varepsilon + (\Delta t)^2). $$
It arises from the triangle inequality and four distinct inequalities,

$$ \| (\mathbf{x}^n, e^n) - (\mathbf{y}^n, g^n) \| \le C\varepsilon, \quad \| (\mathbf{x}, \frac{1}{2}\|\mathbf{v}\|^2) (n\Delta t) - (\mathbf{y}, g)(n\Delta t) \| \le C\varepsilon, $$

$$ \varepsilon \|\mathbf{w}^n - \mathbf{v}(n\Delta t)\| \le C\varepsilon, \quad \| (\mathbf{y}^n, g^n) - (\mathbf{y}, g)(n\Delta t) \| \le C(\Delta t)^2. $$

The first inequality stems from Proposition 4.3, the second one from Proposition 2.3, the third one from uniform bounds on $(\mathbf{w}^n)$ and $\mathbf{v}$ (the former being part of Proposition 3.2) and the fourth is a direct convergence analysis for the Scheme (1.4) towards System (1.2) that we omit as it is a non-stiff simpler version of the arguments detailed to prove Theorem 5.7.

## APPENDIX A. IMPLICIT SOLVING

The asymptotic and convergence analyses of (1.4) are carried out independently of knowing whether the implicit schemes (1.4) and (1.5) do possess unique solutions. Our present goal is to briefly check that this may be enforced by taking $\Delta t$ small (independently of $\varepsilon$).

**Proposition A.1.** Let $T > 0$, $M > 0$. There exist $\delta_0 > 0$ such that when $0 < \Delta t \le \delta_0$ if $(\mathbf{y}, g)$ satisfies $|g| \le \frac{1}{2}M^2$ then there exists a unique $(\mathbf{y}', g')$ such that

$$ (A.1) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\mathbf{y}' - \mathbf{y}}{\Delta t} = -\frac{\mathbf{E}^{\perp}}{b} \left( \frac{\mathbf{y}' + \mathbf{y}}{2} \right) - \frac{g' + g}{2} \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \left( \frac{\mathbf{y}' + \mathbf{y}}{2} \right), \\ \\ \displaystyle \frac{g' - g}{\Delta t} = -\frac{\phi(\mathbf{y}') - \phi(\mathbf{y})}{\Delta t}. \end{array} \right. $$

*Proof.* The second equation provides $g'$ as a function of $(\mathbf{y}, g, \mathbf{y}')$ that may be inserted in the first equation so as to reduce solving System A.1 to finding a solution $\mathbf{y}'$ to a single scalar equation, explicitly written as

$$ \mathbf{y}' = \mathbf{y} - \Delta t \left( \frac{\mathbf{E}^{\perp}}{b} \left( \frac{\mathbf{y}' + \mathbf{y}}{2} \right) - \left( g - \frac{\phi(\mathbf{y}') - \phi(\mathbf{y})}{2} \right) \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \left( \frac{\mathbf{y}' + \mathbf{y}}{2} \right) \right). $$

With respect to the $\mathbf{y}'$ variable the right-hand side of the foregoing equation has Lipschitz constant $O(\Delta t)$ (when $g$ is constrained to vary in a compact set). Thus the proposition follows from the Banach fixed-point theorem. $\square$

**Proposition A.2.** Let $T > 0$, $M > 0$, $\varepsilon_0 > 0$. There exist $\delta_0 > 0$ such that when $0 < \varepsilon \le \varepsilon_0$ and $0 < \Delta t \le \delta_0$ if $(\mathbf{x}, e, \mathbf{w})$ satisfies $\|\mathbf{w}\| \le M$, $|e| \le \frac{1}{2}M^2$ then there exists a unique $(\mathbf{x}', e', \mathbf{w}')$ such that

$$ (A.2) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\mathbf{x}' - \mathbf{x}}{\Delta t} = \frac{\mathbf{w}' + \mathbf{w}}{2\varepsilon} - \left( \frac{e' + e}{2} - \frac{1}{2} \left\| \frac{\mathbf{w}' + \mathbf{w}}{2} \right\|^2 \right) \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right), \\ \\ \displaystyle \frac{e' - e}{\Delta t} = -\frac{\phi(\mathbf{x}') - \phi(\mathbf{x})}{\Delta t}, \\ \\ \displaystyle \frac{\mathbf{w}' - \mathbf{w}}{\Delta t} = \frac{1}{\varepsilon} \mathbf{E} \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) - b \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \frac{(\mathbf{w}')^{\perp} + \mathbf{w}^{\perp}}{2\varepsilon^2}. \end{array} \right. $$

*Proof.* The strategy is to make the most of stability properties in Lemma 3.1 so as to extend the proof of Proposition A.1.

To begin with we observe that one may replace $\mathbf{w}'$ with $\mathbf{w} := (\mathbf{w}'+\mathbf{w})/2$ since $\mathbf{w}'$ is then recovered through $\mathbf{w}' = \mathbf{w} + 2(\mathbf{w}-\mathbf{w})$. Before repeating the elements of the proof of Proposition A.1 we need to derive an a priori bound on $\mathbf{w}$. With $\mathbf{J}$ as in (3.1) (and Lemma 3.1) and $\lambda := \Delta t/\varepsilon^2$ the third equation of (A.2) may be written as

$$ (A.3) \qquad \mathbf{w} = \left( \mathrm{Id} + \frac{\lambda}{2} b \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \mathbf{J} \right)^{-1} \left( \mathbf{w} + \frac{\varepsilon \lambda}{2} \mathbf{E} \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \right) $$
so that it follows from Lemma 3.1 that for some constant $K$ (depending on the lower bound on $b$) any solution $(\mathbf{x}', e', \mathbf{w}')$ to (A.2) satisfies $\|\mathbf{w}\| \le K (M + \varepsilon_0 \|\mathbf{E}\|_{L^\infty})$.

Then we note that the second equation of System A.2 provides $e'$ as a function of $(\mathbf{x}, e, \mathbf{x}')$ and that the outcome may be inserted in its first equation so as to reduce it to

$$ \mathbf{x}' = \mathbf{x} + \Delta t \left( \frac{\mathbf{w}}{\varepsilon} - \left( e - \frac{\phi(\mathbf{x}') - \phi(\mathbf{x})}{2} - \frac{1}{2} \|\mathbf{w}\|^2 \right) \nabla_{\mathbf{x}}^{\perp} \left( \frac{1}{b} \right) \left( \frac{\mathbf{x}' + \mathbf{x}}{2} \right) \right). $$

Heeding the a priori bound on $\mathbf{w}$, we observe that with respect to the $\mathbf{x}'$ variable the right-hand side of the foregoing equation has Lipschitz constant $O(\Delta t)$. Thus provided that $\Delta t$ is sufficiently small (independently of $\varepsilon$) the Banach fixed-point theorem implies that this equation prescribes uniquely $\mathbf{x}'$ as a function of $(\mathbf{x}, e, \mathbf{w})$. For latter use we note that with respect to the $\mathbf{w}$ variable the Lipschitz constant of the implicit map $(\mathbf{x}, e, \mathbf{w}) \mapsto \mathbf{x}'$ is of size $O(\Delta t/\varepsilon)$.

System (A.2) is thus reduced to Equation (A.3) where $\mathbf{x}'$ is replaced with the foregoing function of $(\mathbf{x}, e, \mathbf{w})$. With this point of view Lemma 3.1 implies that, with respect to the $\mathbf{z}$ variable the right-hand side of Equation (A.3) has Lipschitz constant $O(\Delta t)$. Therefore applying the Banach fixed-point theorem again concludes the proof. $\square$

Note that Propositions A.1 and A.2 deal with the solving of a single step of respective schemes. Nevertheless they may be combined with uniform estimates, as provided by Proposition 3.2, to be turned into verifications of the solvability property until times $0 \le n\Delta t \le T$ provided that $\Delta t$ is sufficiently small (depending on $T$ but not on $\varepsilon$).

## 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] P. Chartier, N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field. *Math. Comp.*, 88(320):2697-2736, 2019.

[4] P. Degond and F. Deluzet. Asymptotic-preserving methods and multiscale models for plasma physics. *J. Comput. Phys.*, 336:429-457, 2017.

[5] F. Filbet and L. M. Rodrigues. Asymptotically stable particle-in-cell methods for the Vlasov-Poisson system with a strong external magnetic field. *SIAM J. Numer. Anal.*, 54(2):1120-1146, 2016.

[6] F. Filbet and L. M. Rodrigues. Asymptotically preserving particle-in-cell methods for inhomogeneous strongly magnetized plasmas. *SIAM J. Numer. Anal.*, 55(5):2416-2443, 2017.

[7] F. Filbet and L. M. Rodrigues. Asymptotics of the three-dimensional Vlasov equation in the large magnetic field limit. *J. Éc. polytech. Math.*, 7:1009-1067, 2020.

[8] F. Filbet and L. M. Rodrigues. Asymptotically preserving particle methods for strongly magnetized plasmas in a torus. *J. Comput. Phys.*, 480:Paper No. 112015, 23, 2023.

[9] F. Filbet, L. M. Rodrigues, and K. H. Trinh. A modified Crank-Nicolson scheme for the Vlasov-Poisson system with a strong external magnetic field. *arXiv:2507.02459*, 2025.

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

[11] E. Hairer, C. Lubich, and Y. Shi. Large-stepsize integrators for charged-particle dynamics over multiple time scales. *Numer. Math.*, 151(3):659-691, 2022.

[12] E. Hairer, C. Lubich, and B. Wang. A filtered Boris algorithm for charged-particle dynamics in a strong magnetic field. *Numer. Math.*, 144(4):787-809, 2020.

[13] E. Hairer, C. Lubich, and G. Wanner. *Geometric numerical integration*, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. Structure-preserving algorithms for ordinary differential equations.

[14] E. Hairer, S. P. Nørsett, and G. Wanner. *Solving ordinary differential equations. I*, volume 8 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 1993. Nonstiff problems.

[15] E. Hairer and G. Wanner. *Solving ordinary differential equations. II*, volume 14 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 1996. Stiff and differential-algebraic problems.

[16] S. Jin. Asymptotic-preserving schemes for multiscale physical problems. *Acta Numer.*, 31:415-489, 2022.
[17] K. Miyamoto. *Plasma physics and controlled nuclear fusion*, volume 38 of Springer Series on Atomic, Optical, and Plasma Physics. Springer-Verlag Berlin-Heidelberg, 2006.

[18] B. Wang and Y. Jiang. Semi-discretization and full-discretization with improved accuracy for charged-particle dynamics in a strong nonuniform magnetic field. *ESAIM Math. Model. Numer. Anal.*, 57(4):2427-2450, 2023.

[19] B. Wang and X. Zhao. Error estimates of some splitting schemes for charged-particle dynamics under strong magnetic field. *SIAM J. Numer. Anal.*, 59(4):2075-2105, 2021.

[20] R. Yin. Convergence of some conservative Strang splitting methods for charged-particle dynamics under a strong magnetic field. *J. Comput. Appl. Math.*, 460:Paper No. 116430, 18, 2025.

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

UNIVERSITÉ DE RENNES, INSTITUT DE RECHERCHE EN MATHÉMATIQUES DE RENNES, F-35000 RENNES, FRANCE
*Email address: luis-miguel.rodrigues@univ-rennes.fr*

UNIVERSITÉ DE RENNES, INSTITUT DE RECHERCHE EN MATHÉMATIQUES DE RENNES, F-35000 RENNES, FRANCE
*Email address: kim-han.trinh@univ-rennes.fr*