# Uniformly higher order accurate schemes for dynamics of charged particles under fast oscillating magnetic fields

Megala Anandan, Benjamin Boutin, Nicolas Crouseilles

► To cite this version:

Megala Anandan, Benjamin Boutin, Nicolas Crouseilles. Uniformly higher order accurate schemes for dynamics of charged particles under fast oscillating magnetic fields. IMA Journal of Numerical Analysis, 2025, 048, pp.1-26. hal-04773938

HAL Id: hal-04773938

https://inria.hal.science/hal-04773938v1

Submitted on 8 Nov 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Uniformly higher order accurate schemes for dynamics of charged
particles under fast oscillating magnetic fields*

Megala Anandan.¹, Benjamin Boutin², and Nicolas Crouseilles³

¹Indian Institute of Science, C.V. Raman Road, 560012, Bangalore, India. mail: megalaa@iisc.ac.in.

²Univ Rennes, CNRS, IRMAR - UMR CNRS 6625, F-35000 Rennes, France. mail: benjamin.boutin@univ-rennes.fr.

³Université de Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042 Rennes, France & ENS Rennes. mail: nicolas.crouseilles@inria.fr.

Abstract

This work deals with the numerical approximation of plasmas which are confined by the effect of a
fast oscillating magnetic field (see [12]) in the Vlasov model. The presence of this magnetic field induces
oscillations (in time) to the solution of the characteristic equations. Due to its multiscale character,
a standard time discretization would lead to an inefficient solver. In this work, time integrators are
derived and analyzed for a class of highly oscillatory differential systems. We prove the uniform accuracy
property of these time integrators, meaning that the accuracy does not depend on the small parameter
ε. Moreover, we construct an extension of the scheme which degenerates towards an energy preserving
numerical scheme for the averaged model, when ε → 0. Several numerical results illustrate the capabilities
of the method.

# 1 Introduction

The confinement of charged particles system is of important interest nowadays, in particular due to the
application to fusion devices. The main strategy used in the machines like tokamaks or stellarators is to
apply a strong magnetic field **B**(t, x) (t denotes the time variable and x the space variable) to confine the
particles far from the boundary wall of the device. Indeed, due to the extreme temperature of plasma, it is
crucial to avoid direct contact of the particles with wall material while allowing the plasma to reach fusion
conditions. To confine particles, it is also possible to consider fast oscillating (in time) magnetic fields instead
of strong ones. Indeed, the choice proposed in [12] is the following: **B**(t, x) = θ(t/ε)B(x)b(x), 0 < ε ≪ 1,
where θ is a given P-periodic function, **B**(x) is a scalar positive function, and b(x) is a unitary vector. This
configuration may require less energy to design compared to the strong magnetic field configuration and it
has been proven that this configuration has good confinement properties (see [12]).

In this work, we are interested in the numerical approximation of the model obtained with this highly-
oscillatory magnetic field. From the Vlasov equation (4) in a four-dimensional setting in phase space, the
external magnetic field is supposed to be **B**(t, x) = θ(t/ε)(0, 0, **B**)$^T$ with **B** > 0. The electric field has to be
decomposed into a gradient and rotational part so that we can formulate the corresponding characteristic
equations as follows:

$$
\dot{x}(t) = v(t), \tag{1}
$$

$$
\dot{v}(t) = E(t, x(t)) + \frac{B}{2\epsilon} \theta'(t/\epsilon) Jx(t) + \theta(t/\epsilon) BJv(t), \quad t \in [0, T] \quad (2)
$$

where $x(t) \in \mathbb{R}^2$ denotes the position, $v(t) \in \mathbb{R}^2$ denotes the velocity of the particle, $E : \mathbb{R}_+ \times \mathbb{R}^2 \to \mathbb{R}^2$ denotes a given irrotational electric field (deriving from an electric potential $E = -\nabla\phi$, with $\phi : \mathbb{R}_+ \times \mathbb{R}^2 \to \mathbb{R}$) and $J$ denotes the canonical $2 \times 2$ symplectic matrix.

Our goal is to construct and analyse efficient numerical schemes to approximate differential systems of
the form (1)-(2). It is worth noticing that the system (1)-(2) is highly oscillatory in time ; in particular the

<footnote>*Funding: The first author was funded by the Ministry of Education, Government of India under the Prime Minister's Research Fellowship/grant PM/MHRD-19-17567.03. The last author benefits from fundings by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 EUROFusion). Views and opinions expressed are however those of the author only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.</footnote>
second order derivatives of $v$ is stiff which makes the derivation of higher order schemes challenging since the consistency error usually involves the time derivatives of the solution. To capture the correct dynamics with standard numerical strategies, the time step $\Delta t$ has to be smaller than $\varepsilon \in (0, 1]$, leading to very expensive methods. Thus, the first goal is to overcome the stiffness induced by the highly oscillatory character of the solution. In addition, the method has to be consistent with the asymptotic (or averaged) model obtained in the limit $\varepsilon \to 0$ (see [12]). These properties are usually referred in the literature as asymptotic-preserving or uniformly accurate strategies which ensures that the limit of the scheme is a scheme for the limit model.

During the last decade, several works have been proposed to solve the highly-oscillatory equations efficiently and accurately. These schemes have been referred as uniformly accurate since the error can be proved to be uniform with respect to $\varepsilon$. One can quote the Duhamel based strategies applied for Klein-Gordon or Dirac equations [9, 14, 13], two-scale methods used for Vlasov equations [22, 18, 25, 24, 26, 17, 23, 16, 30, 15, 21, 33], micro-macro approach [19, 23, 20], integrator scheme for advection-diffusion equations [3, 2], multiscale decomposition or exponential wave integrators for Klein-Gordon, Dirac or Schrodinger equations [8, 4, 6, 7, 5]. In this work, we also consider a Duhamel based approach to design uniformly accurate numerical scheme for (1)-(2). A first order uniformly accurate scheme can be easily obtained by freezing the unknown in time and integrating the highly oscillating part exactly, which mainly consists of integrating the function $\theta$. To derive higher order numerical methods, it is possible to recursively substitute adequately the Duhamel formula in itself, in the spirit of [14, 3].

In the linear context, we will consider the following general framework

$$ \dot{U}(t) = A(t/\varepsilon)U(t) \qquad (3) $$

with $U(t) \in \mathbb{R}^d$ and $A$ a given $d \times d$ matrix-valued $P$-periodic function. In this case, it is possible to derive arbitrarily higher order uniformly accurate explicit schemes. On top of that, a suitable use of the Duhamel formula enables to design (second order) midpoint uniformly accurate numerical schemes. Such schemes are particularly relevant for the configuration where the system (3) has a specific structure (Hamiltonian or A skew-symmetric). Indeed, in this case, the averaged model derived from (3) becomes $\dot{U}(t) = \langle A \rangle \bar{U}(t)$ when $\varepsilon \to 0$, and midpoint method is the method of choice since it preserves some invariants ($L^2$ norm or energy) exactly. Hence, our goal is to propose numerical schemes that are not only uniformly accurate, but also degenerate towards a midpoint scheme for the averaged model when $\varepsilon \to 0$.

Then, we investigate how to generalize the strategy to nonlinear systems of the form

$$ \dot{U}(t) = A(t/\varepsilon)U(t) + g(t, U(t)). $$

Even if the first order case can be easily obtained, the extension to higher order requires more attention. Indeed, substituting the Duhamel formula in the nonlinear term leads to a highly oscillatory nonlinear term in $g$ which is not tractable from a numerical point of view. A suitable expansion allows the oscillatory term to be explicitly revealed, enabling the development of higher order, uniformly accurate explicit schemes in the nonlinear context. However, the extension to midpoint strategy (which usually offers invariant preserving properties) turns out to be too costly in view of Vlasov applications since it involves nonlinear fixed point techniques. To overcome this drawback, we employ the SAV technique introduced in [31] to construct efficient structure-preserving numerical methods for nonlinear problems. After reformulating the system (1)-(2) using the SAV framework, we propose an efficient uniformly accurate midpoint numerical scheme which degenerates (when $\varepsilon \to 0$) towards a midpoint scheme for the averaged model of (1)-(2). It turns out that the averaged model enjoys a Hamiltonian structure (see [12]). Hence, a midpoint (with an explicit complexity) numerical scheme is particularly relevant since it exactly preserves the discrete energy.

The rest of the paper is organized as follows. In Section 2, the models we are interested in are presented, together with their asymptotic limit. In section 3, the numerical schemes are presented in a linear context: higher order explicit schemes and second order midpoint schemes are introduced and their uniform accuracy is proved. In section 4.2, we extend the numerical scheme to the nonlinear case. In Section 5, we explain how these solvers can be integrated in Particle-In-Cell framework to approximate the nonlinear Vlasov-Poisson system. Finally, several numerical results are shown to illustrate the properties of the schemes in different contexts.

## 2 Models

The main motivation comes from a physical application introduced in [12] where the author studies the effect of a fast oscillating external magnetic field on a system of charged particles. Let us thus consider the Vlasov equation with the fast oscillating magnetic field $\mathbf{B}(t) = \theta(t/\varepsilon)(0, 0, B)^T$ (with $B > 0$). A more general model
would include gradient and curvature effects, which will be studied in a future work. But since the magnetic field depends on time, the electric field induced from the Maxwell equations takes the form $\frac{B}{2\epsilon}\theta'(t/\epsilon)Jx$ (see [12]). Hence, the Vlasov equation becomes

$$ \partial_t f^\epsilon + v \cdot \nabla_x f^\epsilon + \left( E(t,x) + \frac{B}{2\epsilon} \theta'(t/\epsilon) Jx + B\theta(t/\epsilon) Jv \right) \cdot \nabla_v f^\epsilon = 0, \quad (4) $$

$$ f^{\epsilon}(0, x, v) = f^{\text{in}}(x, v), \quad (x, v) \in \mathbb{R}^2 \times \mathbb{R}^2, \quad t \in [0, T], \quad (5) $$

where $f^{\epsilon}(t, x, v) \ge 0$ is the distribution function of particles in the phase space and $E = -\nabla\phi$ is a given irrotational electric field. Here, $B$ is a constant, $\theta(s)$ is a given $P$-periodic function of class $C^1$, and $0 < \epsilon \ll 1$. Further, for any $a = (a_1, a_2)^T \in \mathbb{R}^2$, $Ja = J[a_1, a_2]^T = [a_2, -a_1]^T$.

As a transport equation, one has $\frac{d}{dt}f(t, x(t), v(t)) = 0$ which translates the fact $f$ is constant along the characteristic curves. Considering $B=1$, $(x(t), v(t))$ is solution to the following ODE system:

$$ \dot{x} = v, $$

$$ \dot{v} = E(t, x) + \frac{1}{2\epsilon}\theta'(t/\epsilon)Jx + \theta(t/\epsilon)Jv. \quad (6) $$

This system is not only stiff but also involves highly oscillatory terms making its numerical approximation quite challenging since a restrictive condition on the time step is commonly required. Before discussing the construction of appropriate numerical schemes for system (6), let us study its asymptotic behavior. As suggested in [12], we first perform the following change of variable $q := v - \frac{1}{2}\theta(t/\epsilon)Jx$, so as to reformulate (6) as a system in the unknown $(x, q)$:

$$ \dot{x} = q + \frac{1}{2}\theta(t/\epsilon)Jx, \quad (7) $$

$$ \dot{q} = E(t, x) + \frac{1}{2}\theta(t/\epsilon)Jq + \frac{1}{4}\theta^2(t/\epsilon)J^2x. $$

Thanks to this change of variable, the system has the form of a standard highly-oscillatory problem for which the limit $\epsilon \to 0$ is well known [1, 11]. Namely, introducing the averaged value $\langle\theta\rangle$ of $\theta$ (with $\langle\theta\rangle := \frac{1}{P}\int_0^P \theta(s)ds$), the solution $(x, q)$ of (7) converges when $\epsilon \to 0$ towards the solution $(\bar{x}, \bar{q})$ of the limiting system:

$$ \dot{\bar{x}} = \bar{q} + \frac{1}{2}\langle\theta\rangle J\bar{x}, \quad (8) $$

$$ \dot{\bar{q}} = E(t, \bar{x}) + \frac{1}{2}\langle\theta\rangle J\bar{q} + \frac{1}{4}\langle\theta^2\rangle J^2\bar{x}. $$

As observed in [12], this averaged model (8) enjoys a Hamiltonian structure where the Hamiltonian is $H(x, q) = \frac{1}{2}|q|^2 + \frac{1}{2}\langle\theta\rangle q \cdot Jx + \frac{1}{8}\langle\theta^2\rangle|x|^2 + \phi(x)$ so that the system (8) can be reformulated as

$$ \left(\begin{array}{c} \dot{\bar{x}} \\ \dot{\bar{q}} \end{array}\right) = \left(\begin{array}{cc} 0 & I \\ -I & 0 \end{array}\right) \left(\begin{array}{c} \nabla_x H(\bar{x}, \bar{q}) \\ \nabla_q H(\bar{x}, \bar{q}) \end{array}\right), \quad I = \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right). $$

In the rest of the paper, for the sake of simplicity, we will consider the following normalized model which shares the same properties as the 'physical' one (7)

$$ \dot{x} = q + \theta(t/\epsilon)Jx, \quad (9) $$

$$ \dot{q} = E(t, x) + \theta(t/\epsilon)Jq + \theta^2(t/\epsilon)J^2x. $$

For the model (9), the averaged model is characterized by the following (normalized) Hamiltonian which writes

$$ H(\bar{x}, \bar{q}) = \frac{1}{2}|\bar{q}|^2 + \langle\theta\rangle \bar{q} \cdot J\bar{x} + \frac{1}{2}\langle\theta^2\rangle |\bar{x}|^2 + \phi(\bar{x}). \quad (10) $$

In this work, we are interested in the construction of efficient numerical schemes for more general highly oscillatory models of the form

$$ \dot{U} = A(t/\epsilon)U + g(U), \quad U(t=0) = U_0 \in \mathbb{R}^d, \quad (11) $$

where $A$ is $d \times d$ matrix-valued $P$-periodic function and $g: \mathbb{R}^d \to \mathbb{R}^d$ is a smooth given function. It is known (see [1, 11]) that the solution $U(t)$ of (11) converges towards $\bar{U}(t)$ which is the solution of

$$ \dot{\bar{U}} = \langle A \rangle \bar{U} + g(\bar{U}), \quad \text{where } \langle A \rangle = \frac{1}{P} \int_0^P A(s)ds \quad \text{and} \quad \bar{U}(t=0) = \bar{U}_0. \quad (12) $$
As an example, the system (9) involves the following oscillating 4 × 4 matrix valued function A(t/ε) and its average 〈 A 〉 defined by

$$
A(t/\epsilon) = \begin{pmatrix} \theta(t/\epsilon)J & I \\ -\theta^2(t/\epsilon)I & \theta(t/\epsilon)J \end{pmatrix} \quad \text{and} \quad \langle A \rangle = \begin{pmatrix} \langle \theta \rangle J & I \\ -\langle \theta^2 \rangle I & \langle \theta \rangle J \end{pmatrix}. \tag{13}
$$

Let us remark that in the linear case g = 0, the averaged model (12) reduces to a simple linear system

$\dot{\bar{U}} = \langle A \rangle \bar{U}$, whose dynamics is determined by the eigenvalues of $\langle A \rangle$. This will be useful to compare with the numerical results.

3 Uniformly accurate numerical schemes: linear case

In this section, we present numerical discretizations of highly oscillatory linear ODEs of the form

$$
\dot{U}(t) = A(t/\epsilon)U(t), \quad U(t=0) = U_0, \tag{14}
$$

where $U : \mathbb{R}_+ \mapsto \mathbb{R}^d$, $U_0 \in \mathbb{R}^d$ and $A : \mathbb{R}_+ \mapsto \mathcal{M}_{d,d}(\mathbb{R})$ is $P$-periodic. Let us recall that when $\epsilon$ tends to 0, the solution to (14) converges towards the solution to the following averaged model

$$
\dot{\bar{U}}(t) = \langle A \rangle \bar{U}(t), \quad \bar{U}(t=0) = U_0. \tag{15}
$$

For the time discretization, we will consider the time step $\Delta t > 0$ and the following notations for $n \in \mathbb{N}$:
$t_n = n\Delta t$ and $U_n \in \mathbb{R}^d$ that denotes the numerical solution at time $t_n$, hopefully an approximation of $U(t_n)$.

3.1 Arbitrarily higher order accurate scheme: explicit case

Before presenting the general framework, let us introduce the Duhamel formula for the solution of (14) step
by step, by integrating it between t_n and t_{n+1}:

$$
U(t_{n+1}) = U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) U(s_1) ds_1. \quad (16)
$$

Integrating (14) again, but between t_n and s_1 gives:

$$
U(s_1) = U(t_n) + \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2,
$$

and substituting this in (16) leads to

$$
\begin{align*}
U(t_{n+1}) &= U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) [U(t_n) + \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2] ds_1 \\
&= U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) ds_1 U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2 ds_1.
\end{align*}
$$

Let $p \ge 1$ be a fixed integer, it is of course possible to recursively substitute $(p-1)$ times the Duhamel formula into itself to get

$$
\begin{align*}
U(t_{n+1}) &= U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) ds_1 U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2 ds_1 \\
&= U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) ds_1 U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) ds_2 ds_1 U(t_n) \\
&\quad + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) \int_{t_n}^{s_2} A(s_3/\epsilon) U(s_3) ds_3 ds_2 ds_1 \\
&= U(t_n) + \sum_{k=1}^{p-1} \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{k-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_k/\epsilon) ds_k ds_{k-1} \cdots ds_1 U(t_n) \\
&\quad + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{p-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_p/\epsilon) U(s_p) ds_p ds_{p-1} \cdots ds_1. \tag{17}
\end{align*}
$$
Thus, the nested Duhamed formula (17) suggests the following p-order numerical scheme for $U_n \approx U(t_n)$:

$$ U_{n+1} = U_n + \sum_{k=1}^{p-1} \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{k-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_k/\epsilon) ds_k ds_{k-1} \cdots ds_1 U_n \\ + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{p-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_p/\epsilon) ds_p ds_{p-1} \cdots ds_1 U_n, \quad (18) $$

$$ = U_n + \left( \sum_{k=1}^{p} H_k \right) U_n, \quad \text{with} \quad H_k = \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{k-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_k/\epsilon) ds_k ds_{k-1} \cdots ds_1. \quad (19) $$

**Remark 1.** When required, the integrals (with respect to $t$) of products of the periodic function $A(t/\epsilon)$, that are involved in the matrices $H_m$, are computed analytically. When difficult or even impossible, the function can be expressed through Fourier series $(\theta(t/\epsilon))^m = \sum_{k \in \mathbb{Z}} \hat{\theta}_{mk} e^{ikt/\epsilon}$ and suitable quadrature $(\theta(t/\epsilon))^m \simeq \sum_{|k| \le K} \hat{\theta}_{mk} e^{ikt/\epsilon}$ can be used for their approximations.

The numerical scheme (18) is uniformly accurate of order $p$, as stated in the following theorem.

**Theorem 1.** Let $U(t=0) \in \mathbb{R}^d$. and $A(s=t/\epsilon) \in \mathcal{M}_{d,d}(\mathbb{R})$ be a given bounded P-periodic function of $s$. Consider the solution $t \mapsto U(t) \in \mathbb{R}^d$ to (14) over a given time interval $[0,T]$. The numerical scheme (18) with initial data $U_0 = U(t=0)$ is uniformly accurate of order $p$: for $\Delta t > 0$, denote $N$ the largest integer such that $t_N \le T$, then

$$ \|U(t_n) - U_n\| \le C \Delta t^p $$

for all $n = 1, 2, \dots, N$, with a constant $C$ independent of $\Delta t$ and $\epsilon$.

*Proof.* From Lemma B-1, $U(t)$ is bounded over the time interval $[0, T]$, uniformly in $\epsilon$. We consider the difference between (17) and (18) to obtain an induction on the error, as in [3]

$$ \begin{aligned} e_{n+1} &= e_n + \sum_{m=1}^{p-1} H_m e_n + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{p-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_p/\epsilon) (U(s_p) - U_n) ds_p ds_{p-1} \cdots ds_1 \\ &= \left( I + \sum_{m=1}^{p-1} H_m + H_p \right) e_n \\ &\quad + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_{p-1}} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_p/\epsilon) (U(s_p) - U(t_n)) ds_p ds_{p-1} \cdots ds_1 \\ &= \left( I + \sum_{m=1}^{p} H_m \right) e_n \\ &\quad + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \cdots \int_{t_n}^{s_p} \int_{t_n}^{s_p} A(s_1/\epsilon) A(s_2/\epsilon) \cdots A(s_p/\epsilon) A(s_{p+1}/\epsilon) U(s_{p+1}) ds_{p+1} ds_p ds_{p-1} \cdots ds_1 \end{aligned} $$

where $t_n \le s_{p+1} \le s_p$. Further since $A(t/\epsilon)$ and $U(t)$ are uniformly bounded in $\epsilon$, we may denote the finite quantities $C_m = \|A\|_{L^\infty}^m$ for $1 \le m \le p$, so as to bound the norm of $e_{n+1}$ by $\|e_{n+1}\| \le (1 + \sum_{m=1}^p C_m \Delta t^m) \|e_n\| + C \Delta t^{p+1}$, and we conclude by the Gronwall lemma. $\square$

In addition to the previous theorem, the following proposition identifies the asymptotic numerical scheme, which is obtained by taking the limit in (18). As expected, thanks to the UA property, this asymptotic scheme is a p-th order accurate scheme for (15).

**Proposition 1.** The numerical scheme (18) converges when $\epsilon \to 0$ to the p-th order accurate scheme below, consistent with the limiting equation (15):

$$ U_{n+1} = \left( \sum_{k=0}^{p} \frac{\Delta t^k}{k!} \langle A \rangle^k \right) U_n. $$

*Proof.* Using lemma B-2, we first get $\int_{t_n}^{s_{m-1}} A(t/\epsilon)dt = (s_{m-1} - t_n)\langle A \rangle + \mathcal{O}(\epsilon)$ and then

$$ \begin{aligned} \int_{t_n}^{s_{m-2}} \int_{t_n}^{s_{m-1}} A(s_{m-1}/\epsilon) A(s_m/\epsilon) ds_m ds_{m-1} &= \int_{t_n}^{s_{m-2}} A(s_{m-1}/\epsilon) \int_{t_n}^{s_{m-1}} A(s_m/\epsilon) ds_m ds_{m-1} \\ &= \int_{t_n}^{s_{m-2}} A(s_{m-1}/\epsilon)(s_{m-1} - t_n) ds_{m-1} \langle A \rangle + \mathcal{O}(\epsilon) = \frac{(s_{m-2} - t_n)^2}{2} \langle A \rangle^2 + \mathcal{O}(\epsilon). \end{aligned} $$
By induction, we similarly obtain

$$
\begin{align*}
H_m &= \int_{t_n}^{t_{n+1}} \int_{t_n}^{s_1} \dots \int_{t_n}^{s_{m-1}} A(s_1/\epsilon) A(s_2/\epsilon) \dots A(s_m/\epsilon) ds_m ds_{m-1} \dots ds_1 \\
&= \int_{t_n}^{t_{n+1}} \frac{(s_1 - t_n)^{m-1}}{(m-1)!} ds_1 \langle A \rangle^m + \mathcal{O}(\epsilon) = \frac{\Delta t^m}{m!} \langle A \rangle^m + \mathcal{O}(\epsilon).
\end{align*}
$$

Therefore, the scheme (18) becomes when $\epsilon \to 0$: $U_{n+1} = U_n + (\sum_{k=1}^p \frac{\Delta t^k}{k!} \langle A \rangle^k) U_n + \mathcal{O}(\epsilon)$. Passing to the limit $\epsilon \to 0$ leads to the result. $\square$

3.2 Uniformly accurate numerical schemes: midpoint case

In this part, we investigate midpoint schemes combined with the strategy presented in the previous subsection.
The reason for studying a possible extension of explicit schemes to midpoint schemes lies in the structure of
the asymptotic model. Indeed, under some assumptions on the matrix $\langle A \rangle$ (skew symmetry or hamiltonian
structure), the asymptotic model preserves some invariants ($L^2$ norm or some energies) and midpoint schemes
are known to preserve the discrete counterpart. Hence, our goal is to construct a numerical scheme that is
uniformly accurate and degenerates towards a scheme (for the averaged model) that also preserves the discrete
invariants, as $\epsilon \to 0$.

Before introducing our strategy, let us first discuss a direct (but naive) way to construct midpoint schemes
for (14). A direct extension of the explicit scheme presented in the previous subsection would be

$$
U_{n+1} = U_n + \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \frac{U_{n+1} + U_n}{2}. \quad (20)
$$

However, even if such a scheme converges towards a midpoint scheme for the averaged model (15), it is only
first order uniformly accurate (even if it is second order accurate for any fixed $\epsilon$). See Appendix A.

We now improve the previous scheme (20) so as to get a second order UA scheme that degenerates towards
a midpoint scheme for the averaged model (8). On one side, we have the forward Duhamel formula

$$
U(s) = U(t_n) + \int_{t_n}^{s} A(s_1/\epsilon)U(s_1)ds_1, \text{ for } s \ge t_n, \quad (21)
$$

and on the other side, we have the backward Duhamel formula

$$
U(s) = U(t_{n+1}) - \int_{s}^{t_{n+1}} A(s_{1}/\epsilon)U(s_{1})ds_{1}, \text{ for } s \le t_{n+1}. \quad (22)
$$

Looking back at the derivation of the UA scheme of previous subsection, we write

$$
U(t_{n+1}) = U(t_n) + \int_{t_n}^{t_{n+1}} A(s/\epsilon)U(s)ds, \quad (23)
$$

and the strategy was to replace $U(s)$ by the corresponding Duhamel formula which enabled us to get explicit
UA schemes. To derive midpoint UA schemes, we now instead replace $U(s)$ by the average of (21) and (22):

$$
U(s) = \frac{U(t_{n+1}) + U(t_n)}{2} + \frac{1}{2} \int_{t_n}^{s} A(s_1/\epsilon)U(s_1)ds_1 - \frac{1}{2} \int_{s}^{t_{n+1}} A(s_1/\epsilon)U(s_1)ds_1, \quad (24)
$$

so that (23) becomes

$$
\begin{equation}
\begin{split}
U(t_{n+1}) = U(t_n) &+ \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \frac{U(t_{n+1}) + U(t_n)}{2} \\
& + \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{t_n}^{s} A(s_1/\epsilon) U(s_1) ds_1 ds - \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{s}^{t_{n+1}} A(s_1/\epsilon) U(s_1) ds_1 ds.
\end{split}
\tag{25}
\end{equation}
$$

Considering the numerical unknown $U_n \approx U(t_n)$, one has the following second order UA numerical scheme

$$
\begin{equation}
\begin{split}
U_{n+1} = U_n &+ \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \frac{U_{n+1} + U_n}{2} \\
&+ \frac{1}{2} \left( \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{t_n}^{s} A(s_1/\epsilon) ds_1 ds - \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{s}^{t_{n+1}} A(s_1/\epsilon) ds_1 ds \right) \frac{U_{n+1} + U_n}{2}.
\end{split}
\tag{26}
\end{equation}
$$
**Theorem 2.** Let $U(t=0) \in \mathbb{R}^d$ and $A(s=t/\varepsilon) \in \mathcal{M}_{d,d}(\mathbb{R})$ be a given bounded P-periodic function of s. Consider the solution $t \mapsto U(t) \in \mathbb{R}^d$ to (14) over a given time interval $[0, T]$. The numerical scheme (26) with initial data $U_0 = U(t=0)$ is second order uniformly accurate: for $\Delta t > 0$, denote $N$ the largest integer such that $t_N \le T$, then

$$ \|U(t_n) - U_n\| \le C \Delta t^2 $$

for all $n = 1, 2, \dots, N$, with a constant $C$ independent of $\Delta t$ and of $\varepsilon$.

*Proof.* Here again, from Lemma B-1, $U(t)$ is bounded over the time interval $[0, T]$, uniformly in $\varepsilon$. Let us find an induction formula on the error $e_{n+1} = U(t_{n+1}) - U_{n+1}$ by considering the difference between (25) and (26)

$$ \begin{aligned} e_{n+1} = e_n &+ \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) ds [e_{n+1} + e_n] + \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) [U(s_1) - \frac{U(t_{n+1}) + U(t_n)}{2}] ds_1 ds \\ & - \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{s}^{t_{n+1}} A(s_1/\varepsilon) [U(s_1) - \frac{U(t_{n+1}) + U(t_n)}{2}] ds_1 ds \\ & + \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) [\frac{e_{n+1} + e_n}{2}] ds_1 ds - \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{s}^{t_{n+1}} A(s_1/\varepsilon) [\frac{e_{n+1} + e_n}{2}] ds_1 ds. \end{aligned} $$

From (24), the third term (and similarly for the fourth term) becomes

$$ \begin{aligned} & \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) [U(s_1) - \frac{U(t_{n+1}) + U(t_n)}{2}] ds_1 ds \\ & = \frac{1}{4} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) [\int_{t_n}^{s_1} A(s_2/\varepsilon) U(s_2) ds_2 - \int_{s_1}^{t_{n+1}} A(s_2/\varepsilon) U(s_2) ds_2] ds_1 ds \end{aligned} $$

so that the error induction becomes

$$ \begin{aligned} e_{n+1} = e_n &+ \left( \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) ds + \frac{1}{4} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) ds_1 ds - \frac{1}{4} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{s}^{t_{n+1}} A(s_1/\varepsilon) ds_1 ds \right) (e_{n+1} + e_n) \\ & + \frac{1}{4} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) \left[ \int_{t_n}^{s_1} A(s_2/\varepsilon) U(s_2) ds_2 - \int_{s_1}^{t_{n+1}} A(s_2/\varepsilon) U(s_2) ds_2 \right] ds_1 ds \\ & - \frac{1}{4} \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{s}^{t_{n+1}} A(s_1/\varepsilon) \left[ \int_{t_n}^{s_1} A(s_2/\varepsilon) U(s_2) ds_2 - \int_{s_1}^{t_{n+1}} A(s_2/\varepsilon) U(s_2) ds_2 \right] ds_1 ds. \end{aligned} $$

From the uniform boundedness of $U$, the terms in the two last lines are uniformly bounded by $C\Delta t^3$ since the intervals of integration are all smaller than $\Delta t$. One can then write

$$ \|e_{n+1}\| \le \|e_n\| + (C\Delta t + C\Delta t^2)(\|e_{n+1}\| + \|e_n\|) + C\Delta t^3, $$

or $(1 - C\Delta t) \|e_{n+1}\| \le (1 + C\Delta t) \|e_n\| + C\Delta t^3$. Thus, we obtain (for $\Delta t < 1/C$)

$$ \|e_{n+1}\| \le \frac{1+C\Delta t}{1-C\Delta t} \|e_n\| + \frac{C}{1-C\Delta t} \Delta t^3 \le (1+K\Delta t) \|e_n\| + C(1+K\Delta t) \Delta t^3 \le (1+K\Delta t) \|e_n\| + K\Delta t^3, $$

and we use the Gronwall lemma to conclude that the error is uniformly bounded by $C\Delta t^2$. We used the simple estimates $(1 + C\Delta t)/(1 - C\Delta t) \le 1 + K\Delta t$ and $1/(1 - C\Delta t) \le 1 + K\Delta t$. $\square$

As done in the previous subsection, we now identify the asymptotic numerical scheme, that is the numerical scheme obtained by passing to the limit in (26).

**Proposition 2.** The midpoint scheme (26) converges when $\varepsilon \to 0$ to the following numerical scheme

$$ \bar{U}_{n+1} = \bar{U}_n + \Delta t \langle A \rangle \frac{\bar{U}_n + \bar{U}_{n+1}}{2}, \qquad (27) $$

which is a second order approximation of the solution $\bar{U}(t_{n+1}) = \exp(\langle A \rangle (n+1)\Delta t)U_0$ of (15). Moreover, when $\langle A \rangle = -\langle A \rangle^T$, the scheme (27) preserves the discrete $L^2$ norm $\|\bar{U}_n\| = \|\bar{U}_0\|$ for all $n \ge 0$. Finally, for $\bar{U}_0 - U_0 = O(\varepsilon)$ and any time $T > 0$, there exists a constant $C > 0$ independent of $\varepsilon$ and $\Delta t$ such that provided $0 \le n\Delta t \le T$, one has

$$ \|\bar{U}_n - U_n\| \le C\varepsilon. $$
Proof. By using Lemma B-2, the scheme (26) is interpreted as

$$
\begin{align*}
U_{n+1} &= U_n + \Delta t \langle A \rangle \frac{U_n + U_{n+1}}{2} + \left( \frac{\Delta t^2}{4} \langle A \rangle^2 - \frac{\Delta t^2}{4} \langle A \rangle^2 \right) \frac{U_n + U_{n+1}}{2} + \mathcal{O}(\Delta t \varepsilon) \\
&= U_n + \Delta t \langle A \rangle \frac{U_n + U_{n+1}}{2} + \mathcal{O}(\Delta t \varepsilon),
\end{align*}
$$

thereby proving the first part of the proposition. For the second part, we consider the scalar product of (27) with (Ūn+1 + Ūn) and use the skew-symmetry of ⟨A⟩. The last point of the proof is obtained by comparing the schemes (26) and (27). Setting

$$
C_n(\varepsilon) = \frac{1}{2} \left( \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{t_n}^{s} A(s_1/\varepsilon) ds_1 ds - \int_{t_n}^{t_{n+1}} A(s/\varepsilon) \int_{s}^{t_{n+1}} A(s_1/\varepsilon) ds_1 ds \right),
$$

and observing then from Lemma B-2, that $C_n(\varepsilon) = \mathcal{O}(\Delta t \varepsilon)$, we have

$$
\begin{align*}
\bar{U}_{n+1} - U_{n+1} &= (\bar{U}_n - U_n) + \Delta t \langle A \rangle \bar{U}_{n+1/2} - \int_{t_n}^{t_{n+1}} A(s/\varepsilon) ds U_{n+1/2} - C_n(\varepsilon) U_{n+1/2} \\
&= (\bar{U}_n - U_n) + \Delta t \langle A \rangle (\bar{U}_{n+1/2} - U_{n+1/2}) + (\Delta t \langle A \rangle - \int_{t_n}^{t_{n+1}} A(s/\varepsilon) ds) U_{n+1/2} - C_n(\varepsilon) U_{n+1/2} \\
&= (I - \Delta t \langle A \rangle / 2)^{-1} \left( (I + \Delta t \langle A \rangle / 2)(\bar{U}_n - U_n) + \mathcal{O}(\Delta t \varepsilon) U_{n+1/2} \right).
\end{align*}
$$

From the boundedness of the numerical solution in bounded times, uniform in $\varepsilon$, and Gronwall inequalities,
we then prove the result.
□

**Remark 2.** The scheme (26) enjoys strong links with Magnus methods [29]. Indeed, in our case, the second order Magnus method would write $U_{n+1} = \exp(\Theta)U_n$, with $\Theta = \int_{t_n}^{t_{n+1}} A(t/\varepsilon)dt + \frac{1}{2} \int_{t_n}^{t_{n+1}} \int_{t_n}^{t} [A(t), A(s)]dsdt$ (with the commutator bracket $[A(t), A(s)] = A(t)A(s) - A(s)A(t)$, which vanishes for scalar $A$ ($d=1$)). A second order expansion of this scheme gives $\exp(\Theta) \approx I + \Theta + \frac{1}{2}\Theta^2$ which corresponds a second order expansion of (26). This link motivates the use of Lie group methods [29] in our context due to their geometric properties in more general applications. As an example, the Magnus scheme $U_{n+1} = \exp(\Theta)U_n$ can be proved to be second order UA, and if $A$ is skew-symmetric, it also preserves the $L^2$ norm of the numerical solution. This class of schemes will be studied in a future work.

3.3 Midpoint schemes for the charged particles model: linear case

In this part, we focus on the charged particles model by applying the scheme (26) to the model (7). Then,
we obtain

$$
\begin{pmatrix} x_{n+1} \\ q_{n+1} \end{pmatrix} = \begin{pmatrix} x_n \\ q_n \end{pmatrix} + \mathcal{B} \left( \frac{x_{n+1}+x_n}{2} \right) + \frac{1}{2} \mathcal{A} \left( \frac{x_{n+1}+x_n}{2} \right), \quad (28)
$$

where $\mathcal{B}$ and $\mathcal{A}$ are defined by

$$
\mathcal{B} = \begin{pmatrix} \int_{t_n}^{t_{n+1}} \theta(s/\varepsilon) ds J & \Delta t I \\ \int_{t_n}^{t_{n+1}} \theta^2(s/\varepsilon) ds J^2 & \int_{t_n}^{t_{n+1}} \theta(s/\varepsilon) ds J \end{pmatrix}, \quad \mathcal{A} = \begin{pmatrix} \mathcal{A}_{1,1} & \mathcal{A}_{1,2} \\ \mathcal{A}_{2,1} & \mathcal{A}_{2,2} \end{pmatrix} \in \mathcal{M}_{4,4}(\mathbb{R}), \quad (29)
$$

where the block components $\mathcal{A}_{i,j} \in \mathcal{M}_{2,2}(\mathbb{R})$ are given by

$$
\mathcal{A}_{1,1} = \int_{t_n}^{t_{n+1}} \theta(s/\varepsilon) \left[ -\int_{t_n}^{s} \theta(s_1/\varepsilon) ds_1 + \int_{s}^{t_{n+1}} \theta(s_1/\varepsilon) ds_1 \right] dsI - \int_{t_n}^{t_{n+1}} \left[ \int_{t_n}^{s} \theta^2(s_1/\varepsilon) ds_1 - \int_{s}^{t_{n+1}} \theta^2(s_1/\varepsilon) ds_1 \right] dsI \quad (30)
$$

$$
\mathcal{A}_{1,2} = \int_{t_n}^{t_{n+1}} \theta(s/\varepsilon) ((s-t_n)-(t_{n+1}-s))dsJ + \int_{t_n}^{t_{n+1}} \left[ \int_{t_n}^{s} \theta(s_1/\varepsilon)ds_1 - \int_{s}^{t_{n+1}} \theta(s_1/\varepsilon)ds_1 \right] dsJ \quad (31)
$$

$$
\mathcal{A}_{2,1} = \int_{t_n}^{t_{n+1}} \theta^2(s/\varepsilon) \left[ -\int_{t_n}^{s} \theta(s_1/\varepsilon) ds_1 + \int_{s}^{t_{n+1}} \theta(s_1/\varepsilon) ds_1 \right] dsJ - \int_{t_n}^{t_{n+1}} \theta(s/\varepsilon) \left[ \int_{t_n}^{s} \theta^2(s_1/\varepsilon) ds_1 - \int_{s}^{t_{n+1}} \theta^2(s_1/\varepsilon) ds_1 \right] dsJ \quad (32)
$$

$$
\mathcal{A}_{2,2} = - \int_{t_n}^{t_{n+1}} \theta^2(s/\varepsilon) ((s-t_n)-(t_{n+1}-s))dsI + \int_{t_n}^{t_{n+1}} \theta(s/\varepsilon) \left[ -\int_{t_n}^{s} \theta(s_1/\varepsilon)ds_1 + \int_{s}^{t_{n+1}} \theta(s_1/\varepsilon)ds_1 \right] dsI. \quad (33)
$$
In the following theorem, we identify the asymptotic numerical scheme by taking the limit $\epsilon \to 0$ in (28). The asymptotic model is satisfied by $(\bar{x}, \bar{q})$ solution to

$$ \dot{\bar{x}} = \bar{q} + \langle \theta \rangle J \bar{x}, \quad (34) $$

$$ \dot{\bar{q}} = \langle \theta \rangle J \bar{q} + \langle \theta^2 \rangle J^2 \bar{x}. \quad (35) $$

This system has a Hamiltonian structure with a Hamiltonian given by (10), and this can be written as $H(x, q) = H_1(x, q) + H_2(x, q)$ defined by

$$ H_1(x, q) = \frac{1}{2} |q|^2 + \frac{1}{2} \langle \theta^2 \rangle |x|^2, \quad \text{and} \quad H_2(x, q) = \langle \theta \rangle q \cdot Jx. \quad (36) $$

Moreover, we have $\frac{d}{dt} H_1(x(t), q(t)) := \frac{d}{dt} (\frac{1}{2}|q(t)|^2 + \frac{1}{2}\langle \theta^2 \rangle |x(t)|^2) = \frac{d}{dt} H_2(x(t), q(t)) = \frac{d}{dt} (\langle \theta \rangle q(t) \cdot Jx(t)) = 0$ where $(x, q)$ is solution to (34)-(35).

**Theorem 3.** The scheme (28) is second order UA and degenerates when $\epsilon \to 0$ towards an energy preserving scheme for the averaged model (34)-(35), whose solution $(\bar{x}_n, \bar{q}_n)_n$ satisfies: $H_1(\bar{x}_{n+1}, \bar{q}_{n+1}) = H_1(\bar{x}_n, \bar{q}_n)$ and $H_2(\bar{x}_{n+1}, \bar{q}_{n+1}) = H_2(\bar{x}_n, \bar{q}_n)$ where $H_1$ and $H_2$ are given by (36).

*Proof.* The UA property follows from the previous theorem. In this proof, we focus on the asymptotic behavior of the scheme (28) when $\epsilon \to 0$. We use the estimates from Lemma B-2 several times, with a crude upper bound of the form $O(\epsilon)$ (instead of $O(\epsilon \Delta t)$ or $O(\epsilon^2)$) so as to highlight the formal limiting scheme. When $\epsilon \to 0$, the matrix $\mathcal{B}$ satisfies

$$ \mathcal{B} = \begin{pmatrix} \int_{t_n}^{t_{n+1}} \theta(s/\epsilon) ds J & \Delta t I \\ \int_{t_n}^{t_{n+1}} \theta^2(s/\epsilon) ds J^2 & \int_{t_n}^{t_{n+1}} \theta(s/\epsilon) ds J \end{pmatrix} = \begin{pmatrix} \Delta t \langle \theta \rangle J & \Delta t I \\ \Delta t \langle \theta^2 \rangle J^2 & \Delta t \langle \theta \rangle J \end{pmatrix} + \mathcal{O}(\epsilon). $$

For the matrix $\mathcal{A}$, let us consider each block component. For the first one, we get from Lemma B-2

$$ \begin{aligned} A_{1,1} &= \left( \int_{t_n}^{t_{n+1}} \theta(s/\epsilon) \int_{t_n}^{s} \theta(s_1/\epsilon) ds_1 ds - \int_{t_n}^{t_{n+1}} \theta(s/\epsilon) \int_{s}^{t_{n+1}} \theta(s_1/\epsilon) ds_1 ds \right. \\ & \qquad \left. + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \theta^2(s_1/\epsilon) ds_1 ds - \int_{t_n}^{t_{n+1}} \int_{s}^{t_{n+1}} \theta^2(s_1/\epsilon) ds_1 ds \right) J^2 \\ &= \frac{\Delta t^2}{2} \langle \theta \rangle^2 - \frac{\Delta t^2}{2} \langle \theta \rangle^2 + \frac{\Delta t^2}{2} \langle \theta^2 \rangle - \frac{\Delta t^2}{2} \langle \theta^2 \rangle + \mathcal{O}(\epsilon) = \mathcal{O}(\epsilon). \end{aligned} $$

For the second component, we get from Lemma B-2

$$ \begin{aligned} A_{1,2} &= \left( \int_{t_n}^{t_{n+1}} \theta(s/\epsilon) \left( (s-t_n) - (t_{n+1}-s) \right) ds + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \theta(s_1/\epsilon) ds_1 ds - \int_{t_n}^{t_{n+1}} \int_{s}^{t_{n+1}} \theta(s_1/\epsilon) ds_1 ds \right) J \\ &= \frac{\Delta t^2}{2} \langle \theta \rangle - \frac{\Delta t^2}{2} \langle \theta \rangle + \frac{\Delta t^2}{2} \langle \theta \rangle - \frac{\Delta t^2}{2} \langle \theta \rangle + \mathcal{O}(\epsilon) = \mathcal{O}(\epsilon). \end{aligned} $$

And similarly, we have for the component $A_{2,1} = \mathcal{O}(\epsilon)$ and $A_{2,2} = \mathcal{O}(\epsilon)$. Hence, the numerical scheme obtained when $\epsilon \to 0$ (with $x_n \to \bar{x}_n$ and $q_n \to \bar{q}_n$) becomes

$$ \begin{pmatrix} \bar{x}_{n+1} \\ \bar{q}_{n+1} \end{pmatrix} = \begin{pmatrix} \bar{x}_n \\ \bar{q}_n \end{pmatrix} + \begin{pmatrix} \Delta t \langle \theta \rangle J & \Delta t I \\ \Delta t \langle \theta^2 \rangle J^2 & \Delta t \langle \theta \rangle J \end{pmatrix} \begin{pmatrix} \frac{\bar{x}_{n+1} + \bar{x}_n}{2} \\ \frac{\bar{q}_{n+1} + \bar{q}_n}{2} \end{pmatrix}. \quad (37) $$

It remains to prove that the scheme (37) preserves the energies $H_1$ and $H_2$.

Let us multiply the first equation of (37) by $\bar{x}_{n+1/2} = (\bar{x}_{n+1} + \bar{x}_n)/2$ to get

$$ (|\bar{x}_{n+1}|^2 - |\bar{x}_n|^2)/\Delta t = \bar{x}_{n+1/2} \cdot \bar{q}_{n+1/2}, $$

since $x \cdot Jx = 0$. Further, multiplying the second equation by $\bar{q}_{n+1/2}$ gives

$$ (|\bar{q}_{n+1}|^2 - |\bar{q}_n|^2)/\Delta t = -\langle \theta^2 \rangle \bar{q}_{n+1/2} \cdot \bar{x}_{n+1/2}. $$

Combining the two last equalities leads to $|\bar{q}_{n+1}|^2 + \langle \theta^2 \rangle |\bar{x}_{n+1}|^2 = |\bar{q}_n|^2 + \langle \theta^2 \rangle |\bar{x}_n|^2$ which corresponds to the $H_1$ preservation. Moreover, multiplying the first equation of (37) by $J$ and considering the scalar product by $\bar{q}_{n+1/2}$ enables to get

$$ \bar{q}_{n+1/2} \cdot J(\bar{x}_{n+1} - \bar{x}_n)/\Delta t = \langle \theta \rangle \bar{q}_{n+1/2} \cdot J J \bar{x}_{n+1/2} + \bar{q}_{n+1/2} \cdot J \bar{q}_{n+1/2} = -\langle \theta \rangle \bar{q}_{n+1/2} \cdot \bar{x}_{n+1/2}. $$
Now, multiplying the second equation of (37) by $J\bar{x}_{n+1/2}$ leads to

$$J\bar{x}_{n+1/2} \cdot (\bar{q}_{n+1} - \bar{q}_n)/\Delta t = \langle \theta \rangle J\bar{x}_{n+1/2} \cdot J\bar{q}_{n+1/2} = \langle \theta \rangle \bar{x}_{n+1/2} \cdot \bar{q}_{n+1/2}.$$

Adding the two last equalities leads to

$$\bar{q}_{n+1/2} \cdot J(\bar{x}_{n+1} - \bar{x}_n) + J\bar{x}_{n+1/2} \cdot (\bar{q}_{n+1} - \bar{q}_n) = 0,$$

but the left hand side can be reformulated as

$$\begin{aligned}
\bar{q}_{n+1/2} \cdot J(\bar{x}_{n+1} - \bar{x}_n) + J\bar{x}_{n+1/2} \cdot (\bar{q}_{n+1} - \bar{q}_n) &= (\bar{q}_{n+1} + \bar{q}_n) \cdot J(\bar{x}_{n+1} - \bar{x}_n) + J(\bar{x}_{n+1} + \bar{x}_n) \cdot (\bar{q}_{n+1} - \bar{q}_n) \\
&= \bar{q}_{n+1} \cdot J\bar{x}_{n+1} - \bar{q}_{n+1} \cdot J\bar{x}_n + \bar{q}_n \cdot J\bar{x}_{n+1} - \bar{q}_n \cdot J\bar{x}_n + J\bar{x}_{n+1} \cdot \bar{q}_{n+1} - J\bar{x}_{n+1} \cdot \bar{q}_n + J\bar{x}_n \cdot \bar{q}_{n+1} - J\bar{x}_n \cdot \bar{q}_n \\
&= \bar{q}_{n+1} \cdot J\bar{x}_{n+1} - \bar{q}_n \cdot J\bar{x}_n,
\end{aligned}
\quad (38)$$

which corresponds to the $H_2$ preservation and concludes the proof. □

**Remark 3.** It is actually sufficient to prove that the scheme (28) degenerates towards a midpoint scheme for the averaged model (34)-(35), when $\varepsilon \to 0$. Indeed, since $H_1$ and $H_2$ are quadratic, we can appeal Theorem 2.1, p101 of [28] to ensure that the midpoint method preserves quadratic invariants. However, the detailed proof will be useful for the nonlinear case.

# 4 Nonlinear system for charged particles under fast oscillating magnetic field

In this section, we consider some possible extensions of the numerical schemes introduced in previous section to the nonlinear case.

## 4.1 Explicit case

We first extend the explicit uniformly accurate numerical schemes presented in subsection 3.1 to nonlinear case. In this subsection, we consider nonlinear ODEs of the form

$$\dot{U}(t) = A(t/\varepsilon)U(t) + g(t, U(t)), \quad U(t=0) = U_0. \quad (39)$$

for which we propose first and second order uniformly accurate schemes.

### 4.1.1 First order accurate scheme

Let us consider the general form for nonlinear systems of the form (39) for which the following first order uniformly accurate numerical scheme can be obtained from the previous section

$$U_{n+1} = U_n + \int_{t_n}^{t_{n+1}} A(t/\varepsilon) dt U_n + \Delta t g_n \text{ with } g_n = g(t_n, U_n). \quad (40)$$

We have the following theorem.

**Theorem 4.** Let $U(t=0) \in \mathbb{R}^d$, $A(s=t/\varepsilon) \in \mathcal{M}_{d,d}(\mathbb{R})$ be a given bounded P-periodic function of s, and g be a given bounded Lipschitz continuous function of $(t, U)$ with Lipschitz constant $K$. Consider the solution $t \mapsto U(t) \in \mathbb{R}^d$ to (39) over a given time interval $[0, T]$. The numerical scheme (40) with initial data $U_0 = U(t=0)$ is first order uniformly accurate:
for $\Delta t > 0$, denote $N$ the largest integer such that $t_N \le T$, then

$$\|U(t_n) - U_n\| \le C\Delta t$$

for all $n = 1, 2, \dots, N$, with a constant $C$ independent of $\Delta t$ and of $\varepsilon$.

*Proof.* Integrating (39) between $t_n$ and $t$ gives,

$$U(t) = U(t_n) + \int_{t_n}^{t} A(s/\varepsilon) U(s) ds + \int_{t_n}^{t} g(s, U(s)) ds.$$
Evaluating the above exact integration at $t = t_{n+1}$, we get

$$U(t_{n+1}) = U(t_n) + \int_{t_n}^{t_{n+1}} A(s/\epsilon) U(s) ds + \int_{t_n}^{t_{n+1}} g(s, U(s)) ds. \quad (41)$$

Then with $e_n = U(t_n) - U_n$, the difference between (40) and (41) gives

$$\begin{aligned}
e_{n+1} &= e_n + \int_{t_n}^{t_{n+1}} A(s/\epsilon) (U(s) - U_n) ds + \int_{t_n}^{t_{n+1}} (g(s, U(s)) - g_n) ds \\
&= e_n + \int_{t_n}^{t_{n+1}} A(s/\epsilon) (U(s) - U(t_n)) ds + \int_{t_n}^{t_{n+1}} A(s/\epsilon) e_n ds + \int_{t_n}^{t_{n+1}} (g(s, U(s)) - g_n) ds \\
&= \left( I + \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \right) e_n + \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{t_n}^{s} \left[ A(s_1/\epsilon) U(s_1) + g(s_1, U(s_1)) \right] ds_1 ds \\
&\quad + \int_{t_n}^{t_{n+1}} (g(s, U(s)) - g_n) ds.
\end{aligned}$$

Further since $\|A(t/\epsilon)\|$, $\|U(t)\|$, and $\|g(t, U(t))\|$ are uniformly bounded and we also have $\|g(s, U(s)) - g_n\| \le K \|(s, U(s)) - (t_n, U_n)\| \le C(\Delta t + \|e_n\|)$, where K is the Lipschitz constant. Then the norm of $e_{n+1}$ becomes

$$\|e_{n+1}\| \le (1 + C\Delta t) \|e_n\| + C\Delta t^2 + C\Delta t (\Delta t + \|e_n\|) \le (1 + C\Delta t) \|e_n\| + C\Delta t^2,$$

where the constant C can changed from one line to the next but it is independent of $\Delta t, \epsilon$. Thus we conclude by Gronwall lemma that $\|e_n\| \le C\Delta t$ with C independent of $\Delta t$ and $\epsilon$, for $n = 1, 2, \dots, N$. $\square$

### 4.1.2 Second order accurate scheme

In this section, we present how to extend the above first order scheme to get a second order uniformly accurate scheme for (39). The following second order uniformly accurate scheme (and the corresponding theorem) holds for $g(t, U(t))$ but for the sake of simplicity, we will present only for the case $g(U(t))$.

To derive a second order scheme (39), we follow the strategy introduced in the linear case: we integrate (39) between $t_n$ and $t_{n+1}$ and perform a recursive substitution to get

$$\begin{aligned}
U(t_{n+1}) &= U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) U(s_1) ds_1 + \int_{t_n}^{t_{n+1}} g(U(s_1)) ds_1 \\
&= U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) ds_1 U(t_n) + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2 ds_1 \\
&\quad + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} g(U(s_2)) ds_2 ds_1 + \int_{t_n}^{t_{n+1}} g(U(t_n) + h(s_1)) ds_1,
\end{aligned} \quad (42)$$

with $h(s_1) := \int_{t_n}^{s_1} [A(s_2/\epsilon) U(s_2) + g(U(s_2))] ds_2$. In view of the numerical scheme, the solution is frozen at time $t_n$ in the third and fourth terms whereas we have to discuss the last term. First, assuming $g$ is smooth, we perform a Taylor expansion of $g(U(t_n) + h(s_1)) \approx g(U(t_n)) + (\nabla g)(U(t_n))h(s_1)$ and second, one uses a suitable approximation $\tilde{h}$ of $h$. We then obtain the following second order numerical discretization of (39)

$$\begin{aligned}
U_{n+1} = U_n &+ \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) ds_1 U_n + \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) ds_2 ds_1 U_n \\
&+ \int_{t_n}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} ds_2 ds_1 g_n + \Delta t g_n + \nabla g|_n \int_{t_n}^{t_{n+1}} \tilde{h}(s_1) ds_1
\end{aligned} \quad (43)$$

where $g_n = g(U_n)$ and $\nabla g|_n = (\nabla g)(U_n)$ and $\tilde{h}$ is an approximation of $h$ given by

$$\tilde{h}(s_1) = \int_{t_n}^{s_1} A(s_2/\epsilon) ds_2 U_n + (s_1 - t_n) g_n.$$

**Theorem 5.** Let $U(t = 0) \in \mathbb{R}^d$, $A(s = t/\epsilon) \in \mathcal{M}_{d,d}(\mathbb{R})$ be a given bounded P-periodic function of s, and $g, \nabla g$ be two given bounded Lipschitz continuous functions of U with Lipschitz constant K. Consider the solution $t \mapsto U(t) \in \mathbb{R}^d$ to (39) over a given time interval $[0, T]$. The numerical scheme (43) with initial data
$U_0 = U(t=0)$ is second order uniformly accurate:

for $\Delta t > 0$, denote $N$ the largest integer such that $t_N \le T$, then

$$ \|U(t_n) - U_n\| \le C\Delta t^2 $$

for all $n = 1, 2, \dots, N$, with a constant $C$ independent of $\Delta t$ and of $\varepsilon$.

*Proof.* Considering the notation $e_n = U(t_n) - U_n$ for $n = 1, 2, \dots, N$, subtracting (42) with (43) gives

$$ \begin{aligned} e_{n+1} &= e_n + \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) ds_1 e_n + \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) \int_{t_n}^{s_1} A(s_2/\varepsilon) (U(s_2) - U_n) ds_2 ds_1 \\ & \qquad + \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) \int_{t_n}^{s_1} (g(U(s_2)) - g_n) ds_2 ds_1 + \int_{t_n}^{t_{n+1}} \left( g(U(t_n) + h(s_1)) - g_n - \nabla g|_n \tilde{h}(s_1) \right) ds_1 \\ &= e_n + I + II + III + IV. \end{aligned} $$

Let us first consider the term $II$:

$$ \begin{aligned} II &= \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) \int_{t_n}^{s_1} A(s_2/\varepsilon) (U(s_2) - U_n) ds_2 ds_1 = \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) \int_{t_n}^{s_1} A(s_2/\varepsilon) (U(s_2) - U(t_n) + e_n) ds_2 ds_1 \\ &= \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) \int_{t_n}^{s_1} A(s_2/\varepsilon) \int_{t_n}^{s_2} \left( A(s_3/\varepsilon) U(s_3) + g(U(s_3)) \right) ds_3 ds_2 ds_1 + \int_{t_n}^{t_{n+1}} A(s_1/\varepsilon) \int_{t_n}^{s_1} A(s_2/\varepsilon) ds_2 ds_1 e_n(44) \end{aligned} $$

From the assumptions on $A$ and $g$, we deduce $\|II\| \le C\Delta t^3 + C\Delta t^2 e_n$. Similarly for the term $III$, we have $\|III\| \le C\Delta t^2 e_n$. Considering the term $IV$, we have

$$ \begin{aligned} IV &= \int_{t_n}^{t_{n+1}} \left( g(U(t_n) + h(s_1)) - g_n - \nabla g|_n \tilde{h}(s_1) \right) ds_1 \\ &= \int_{t_n}^{t_{n+1}} \left( g(U(t_n) + h(s_1)) - g(U(t_n)) - (\nabla g)(U(t_n)) \tilde{h}(s_1) \right) ds_1 \\ & \qquad + \int_{t_n}^{t_{n+1}} \left( g(U(t_n)) - g_n + ((\nabla g)(U(t_n)) - \nabla g|_n) \tilde{h}(s_1) \right) ds_1. \end{aligned} $$

From Lipshitz property on $g$ and $\nabla g$ and from the estimate $\|\tilde{h}(s_1)\| \le C\Delta t$, we deduce that the second term on the right hand side is bounded by $C\Delta t e_n$. Upon Taylor expanding $g(U(t_n) + h(s_1))$ and considering the integral remainder, the first term on the right hand side of above equation becomes

$$ \int_{t_n}^{t_{n+1}} (\nabla g)(U(t_n))(h(s_1) - \tilde{h}(s_1))ds_1 + \int_{t_n}^{t_{n+1}} \int_{0}^{1} (1-\tau)(\text{Hessg})(U(t_n) + \tau h(s_1))(h(s_1), h(s_1))d\tau ds_1 $$

where Hessg denotes the Hessian matrix of $g$ with respect to the variable $U$. Further, we have

$$ \|(h - \tilde{h})(t)\| = \left\| \int_{t_n}^{t} A(s/\varepsilon)(U(s) - U_n)ds + \int_{t_n}^{t} [g(U(s)) - g_n]ds \right\|, $$

but $\|g(U(s)) - g_n\| \le K\|U(s) - U_n\| \le K(e_n + C\Delta t)$ hence $|(h - \tilde{h})(t)| \le C\Delta t^2$. Further we also observe the following for $t_n < s_1 < t_{n+1}$, with the notation $K_1 = \|\text{Hessg}(U)\|$

$$ \|\text{Hessg}(U)(h(s_1), h(s_1))\| \le K_1 \|h(s_1)^2\| \le K_1 \Delta t^2, $$

so that we deduce $\|IV\| \le C\Delta t^3 + C\Delta t e_n$. Then the norm of $e_{n+1}$ becomes: $\|e_{n+1}\| \le (1 + C\Delta t + C\Delta t^2) \|e_n\| + C\Delta t^3$ and we conclude from the Gronwall lemma that $\|e_n\| \le C\Delta t^2$. $\square$

**Remark 4.** Up to some technical and lengthy developments, it is possible to generalize the second order scheme presented above to arbitrary order of accuracy.

## 4.2 Midpoint schemes for the nonlinear charged particles model

To handle the nonlinear models through midpoint schemes, fixed point strategies are usually required. To overcome this drawback, we introduce a scalar auxiliary variable following SAV techniques introduced in [31]. More precisely, considering the system (9) and the averaged model (45) which enjoys a Hamiltonian structure with the Hamiltonian (10), our goal is to design uniformly accurate scheme for (9) which degenerates into an energy preserving scheme for the averaged model (45). Let us remark that contrary to the linear case, we do not have preservation of both Hamiltonians $H_1$ and $H_2$ (defined in (36)).
### 4.2.1 SAV formulation of the averaged model

In the nonlinear case, the averaged model can be written as

$$ 
\begin{aligned}
\dot{\bar{x}} &= \bar{q} + \langle \theta \rangle J \bar{x}, \\
\dot{\bar{q}} &= -\nabla \phi(\bar{x}) + \langle \theta \rangle J \bar{q} + \langle \theta^2 \rangle J^2 \bar{x},
\end{aligned}
\qquad (45) $$

and has a Hamiltonian structure where the Hamiltonian is given by (10). Moreover, it is known that, except in some specific case or when using nonlinear solvers (e.g. averaged gradient method), there is no general numerical scheme preserving the Hamiltonian. We instead consider an extended ODE system based on a modification of the Hamiltonian through the SAV technique [31]. Here, using the additional variable $\bar{r}(t) = \exp(\phi(\bar{x}(t)))$, it is possible to reformulate the system (45) into

$$ 
\begin{aligned}
\dot{\bar{x}} &= \bar{q} + \langle \theta \rangle J \bar{x}, \\
\dot{\bar{q}} &= -\bar{b} + \langle \theta \rangle J \bar{q} + \langle \theta^2 \rangle J^2 \bar{x}, \\
\dot{\bar{r}} &= \bar{r} \bar{b} \cdot \dot{\bar{x}},
\end{aligned}
\qquad (46) $$

where we introduced $\bar{b}(t) = \nabla\phi(\bar{x}(t))$. Other choices for the additional variable $\bar{r}$ are possible, the most used being $\bar{r}(t) = \sqrt{\phi(\bar{x}(t))} + C$ with $C$ large enough. This extended system enjoys an energy preservation property as shown in the next proposition.

**Proposition 3.** For the extended system (46), the following quantity is preserved with time

$$ \bar{H}(\bar{x}, \bar{q}, \bar{r}) = \frac{1}{2} |\bar{q}|^2 + \langle \theta \rangle \bar{q} \cdot J \bar{x} + \frac{1}{2} \langle \theta^2 \rangle |\bar{x}|^2 + \log(\bar{r}). \qquad (47) $$

*Proof.* Multiplying the first equation of (46) by $\langle \theta^2 \rangle \bar{x}$, the second by $\bar{q}$ and combining them with the third equation gives

$$ 
\begin{aligned}
\frac{d}{dt} \left( \frac{1}{2} |\bar{q}|^2 + \frac{1}{2} \langle \theta^2 \rangle |\bar{x}|^2 + \log(\bar{r}) \right) &= \langle \theta^2 \rangle \bar{x} \cdot \bar{q} - \bar{b} \cdot \bar{q} + \langle \theta^2 \rangle J^2 \bar{x} \cdot \bar{q} + \bar{b} \cdot \dot{\bar{x}} \\
&= -\bar{b} \cdot \left( \dot{\bar{x}} - \langle \theta \rangle J \bar{x} \right) + \bar{b} \cdot \dot{\bar{x}} = \bar{b} \cdot \langle \theta \rangle J \bar{x}.
\end{aligned}
 $$

The identity $-\bar{b} \cdot J \bar{x} = \frac{d}{dt}(\bar{q} \cdot J \bar{x})$ follows from

$$ \frac{d}{dt}(\bar{q} \cdot J \bar{x}) = \dot{\bar{q}} \cdot J \bar{x} + \bar{q} \cdot J \dot{\bar{x}} = -\bar{b} \cdot J \bar{x} + \langle \theta \rangle J \bar{q} \cdot J \bar{x} + \langle \theta \rangle \bar{q} \cdot J^2 \bar{x} = -\bar{b} \cdot J \bar{x} $$

□

### 4.2.2 Numerical scheme for the SAV formulation of the averaged model

Next, we propose the following numerical scheme for solving the nonlinear averaged model (46):

$$ 
\begin{aligned}
(\bar{x}_{n+1} - \bar{x}_n) &= \Delta t \bar{q}_{n+1/2} + \Delta t \langle \theta \rangle J \bar{x}_{n+1/2}, \\
(\bar{q}_{n+1} - \bar{q}_n) &= -\Delta t \bar{b} + \Delta t \langle \theta \rangle J \bar{q}_{n+1/2} - \Delta t \langle \theta^2 \rangle \bar{x}_{n+1/2}, \\
\log(\bar{r}_{n+1}) - \log(\bar{r}_n) &= \bar{b} \cdot (\bar{x}_{n+1} - \bar{x}_n),
\end{aligned}
\qquad (48) $$

where $x_{n+1/2} = (x_{n+1} + x_n)/2$ and $q_{n+1/2} = (q_{n+1} + q_n)/2$. The approximation for the quantity $\bar{b}$ (which can be viewed as an approximation of $(\int_{t_n}^{t_{n+1}} \nabla\phi(\bar{x}(s))ds)/\Delta t$) will be discussed later on and has no consequence on the following result:

**Proposition 4.** The numerical solution to the scheme (48) preserves the Hamiltonian $\bar{H}$ given by (47):
$$ \bar{H}(\bar{x}_{n+1}, \bar{q}_{n+1}, \bar{r}_{n+1}) = \bar{H}(\bar{x}_n, \bar{q}_n, \bar{r}_n), \quad n \ge 0. $$

*Proof.* The proof follows the lines of the proof of Proposition 3. Multiplying the first equation of (48) by $\langle \theta^2 \rangle \bar{x}_{n+1/2}$, the second equation by $\bar{q}_{n+1/2}$ and combining them with the third equation gives

$$ 
\begin{aligned}
& \frac{1}{2} \langle \theta^2 \rangle (|\bar{x}_{n+1}|^2 - |\bar{x}_n|^2) + \frac{1}{2} (|\bar{q}_{n+1}|^2 - |\bar{q}_n|^2) + \log(\bar{r}_{n+1}) - \log(\bar{r}_n) = \\
& \Delta t \langle \theta^2 \rangle \bar{x}_{n+1/2} \cdot \bar{q}_{n+1/2} - \Delta t \bar{b} \cdot \bar{q}_{n+1/2} - \Delta t \langle \theta^2 \rangle \bar{x}_{n+1/2} \cdot \bar{q}_{n+1/2} + \bar{b} \cdot (\bar{x}_{n+1} - \bar{x}_n) \\
& = -\bar{b} \cdot \left( (\bar{x}_{n+1} - \bar{x}_n) - \Delta t \langle \theta \rangle J \bar{x}_{n+1/2} \right) + \bar{b} \cdot (\bar{x}_{n+1} - \bar{x}_n) = \Delta t \langle \theta \rangle \bar{b} \cdot J \bar{x}_{n+1/2}.
\end{aligned}
\qquad (49) $$
It remains to check $-\Delta t \bar{b} \cdot J \bar{x}_{n+1/2} = \bar{q}_{n+1} \cdot J \bar{x}_{n+1} - \bar{q}_n \cdot J \bar{x}_n$.

Multiplying the second equation of (48) by $J \bar{x}_{n+1/2}$ leads to

$$ (\bar{q}_{n+1} - \bar{q}_n) \cdot J \bar{x}_{n+1/2} = -\Delta t \bar{b} \cdot J \bar{x}_{n+1/2} + \Delta t \langle \theta \rangle J \bar{q}_{n+1/2} \cdot J \bar{x}_{n+1/2} = -\Delta t \bar{b} \cdot J \bar{x}_{n+1/2} + \Delta t \langle \theta \rangle \bar{q}_{n+1/2} \cdot \bar{x}_{n+1/2}. $$

Considering $\bar{q}_{n+1/2} \cdot J$ multiplied by the first equation of (48) leads to

$$ \bar{q}_{n+1/2} \cdot J(\bar{x}_{n+1} - \bar{x}_n) = \Delta t \langle \theta \rangle \bar{q}_{n+1/2} \cdot J^2 \bar{x}_{n+1/2} = -\Delta t \langle \theta \rangle \bar{q}_{n+1/2} \cdot \bar{x}_{n+1/2}. $$

Adding the two last equalities gives

$$ (\bar{q}_{n+1} - \bar{q}_n) \cdot J \bar{x}_{n+1/2} + \bar{q}_{n+1/2} \cdot J(\bar{x}_{n+1} - \bar{x}_n) = -\Delta t \bar{b} \cdot J \bar{x}_{n+1/2}, $$

and using (38), one gets $-\Delta t \bar{b} \cdot J \bar{x}_{n+1/2} = \bar{q}_{n+1} \cdot J \bar{x}_{n+1} - \bar{q}_n \cdot J \bar{x}_n$ which enables to conclude with (49). $\square$

We end this part by discussing the approximation of $\bar{b}$. The usual choice is $\bar{b} \approx \nabla \phi(\bar{x}(t_{n+1/2}))$ so that second order accuracy is reached through extrapolation techniques [31]:

$$ \bar{b} \approx \nabla \phi(\bar{x}(t_{n+1/2})) \approx -\frac{1}{2} \nabla \phi(\bar{x}_{n-1}) + \frac{3}{2} \nabla \phi(\bar{x}_n). \qquad (50) $$

Even if this choice enables to get a linearly implicit scheme, it requires the knowledge of $x_{n-1}$. Here, another technique is proposed based on the above schemes derived from integration in time. Indeed, we start from

$$ \bar{b} = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \nabla \phi(\bar{x}(s)) ds. \qquad (51) $$

But from the integration on $t \in [t_n, s]$ of the first equation of (46), we get

$$ \bar{x}(s) = \bar{x}(t_n) + \bar{h}(s) \quad \text{with} \quad \bar{h}(s) = \int_{t_n}^{s} \bar{q}(s_1) ds_1 + \int_{t_n}^{s} \langle \theta \rangle J \bar{x}(s_1) ds_1. $$

We can now approximate $\bar{h}$ by $\tilde{h} = (s - t_n)\bar{q}(t_n) + (s - t_n)\langle \theta \rangle J \bar{x}(t_n)$ so that a choice for $\bar{b}$ is obtained by inserting this last approximation in a Taylor expansion for $\nabla \phi$ in (51)

$$ \begin{aligned} \bar{b} &= \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \nabla \phi(\bar{x}(s)) ds = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \nabla \phi(\bar{x}(t_n) + \bar{h}(s)) ds \\ &\approx \nabla \phi(\bar{x}(t_n)) + \nabla^2 \phi(\bar{x}(t_n)) \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \tilde{h}(s) ds \\ &\approx \nabla \phi(\bar{x}_n) + \nabla^2 \phi(\bar{x}_n) \frac{\Delta t}{2} (\bar{q}_n + \langle \theta \rangle J \bar{x}_n), \end{aligned} \qquad (52) $$

which can be proved to be a second order (third locally) approximation of $(\int_{t_n}^{t_{n+1}} \nabla \phi(\bar{x}(s)) ds)/\Delta t$ and only requires the knowledge of $\bar{x}_n$ to update the numerical unknown $(\bar{x}_{n+1}, \bar{q}_{n+1}, \bar{r}_{n+1})$.

### 4.2.3 Uniformly accurate and structure preserving SAV-schemes for the charged particle model

The goal of this part is to propose a numerical scheme for the charged particle model (9) which is second order uniformly accurate and which degenerates when $\varepsilon \to 0$ towards the SAV-midpoint scheme (48) (which preserves the energy thanks to Proposition 4). First, we derive the SAV formulation for (9), inspired by the SAV formulation of the averaged model

$$ \begin{aligned} \dot{x} &= q + \theta Jx, \\ \dot{q} &= -b + \theta Jq + \theta^2 J^2 x, \\ \dot{r} &= rb \cdot \dot{x}. \end{aligned} \qquad (53) $$

To derive a second order UA scheme, let us integrate the first two equations of (53) on $t \in [t_n, t_{n+1}]$ to get

$$ x(t_{n+1}) - x(t_n) = \int_{t_n}^{t_{n+1}} q(t) dt + \int_{t_n}^{t_{n+1}} \theta(t/\varepsilon) Jx(t) dt, \qquad (54) $$

$$ q(t_{n+1}) - q(t_n) = - \int_{t_n}^{t_{n+1}} b(t) dt + \int_{t_n}^{t_{n+1}} \theta(t/\varepsilon) Jq(t) dt + \int_{t_n}^{t_{n+1}} \theta^2(t/\varepsilon) J^2 x(t) dt, \qquad (55) $$
but following the calculations performed in (25), we get

$$
\begin{align*}
x(t) &= \frac{1}{2}(x(t_{n+1} + x(t_n)) + \frac{1}{2}\left[\int_{t_n}^{t}(q(s) + \theta(s/\varepsilon)Jx(s))ds - \int_{t}^{t_{n+1}}(q(s) + \theta(s/\varepsilon)Jx(s))ds\right] \\
q(t) &= \frac{1}{2}(q(t_{n+1} + q(t_n)) - \frac{1}{2}\left[\int_{t_n}^{t}b(s)ds - \int_{t}^{t_{n+1}}b(s)ds\right] \\
&\quad + \frac{1}{2}\left[\int_{t_n}^{t}(\theta(s/\varepsilon)Jq(s) + \theta^2(s/\varepsilon)J^2x(s))ds - \int_{t}^{t_{n+1}}(\theta(s/\varepsilon)Jq(s) + \theta^2(s/\varepsilon)J^2x(s))ds\right].
\end{align*}
$$

Inserting these expressions into (54) and (55), and using midpoint approximation of $(x(s), q(s)) \approx (x_{n+1/2}, q_{n+1/2})$
leads to

$$
x_{n+1} - x_n = (\mathcal{B}_{1,1} + \frac{1}{2}\mathcal{A}_{1,1})x_{n+1/2} + (\mathcal{B}_{1,2} + \frac{1}{2}\mathcal{A}_{1,2})q_{n+1/2} - \frac{1}{2}\int_{t_n}^{t_{n+1}}\left[\int_{t_n}^{t}b(s)ds - \int_{t}^{t_{n+1}}b(s)ds\right]dt, (56)
$$

$$
\begin{equation}
\begin{split}
q_{n+1} - q_n = & - \int_{t_n}^{t_{n+1}} b(s) ds + (\mathcal{B}_{2,1} + \frac{1}{2} \mathcal{A}_{2,1}) x_{n+1/2} + (\mathcal{B}_{2,2} + \frac{1}{2} \mathcal{A}_{2,2}) q_{n+1/2} \\
& - \frac{1}{2} \int_{t_n}^{t_{n+1}} \theta(t/\varepsilon) J \left[ \int_{t_n}^{t} b(s) ds - \int_{t}^{t_{n+1}} b(s) ds \right] dt,
\end{split}
\tag{57}
\end{equation}
$$

where the 2x2 block matrices are given by (30), (31), (32), (33) and the 4x4 matrix $\mathcal{B}$ is given by (29). The $b$
terms with double integrals can be approximated at $s = t_n$ (the last term in (56) vanishes) whereas the term
$-\int_{t_n}^{t_{n+1}} b(s)ds$ will be discussed after. We thus consider the following numerical scheme for (48)

$$
\begin{align}
x_{n+1} - x_n &= (\mathcal{B}_{1,1} + \frac{1}{2}\mathcal{A}_{1,1})x_{n+1/2} + (\mathcal{B}_{1,2} + \frac{1}{2}\mathcal{A}_{1,2})q_{n+1/2}, \nonumber \\
q_{n+1} - q_n &= -\int_{t_n}^{t_{n+1}} b(s)ds + (\mathcal{B}_{2,1} + \frac{1}{2}\mathcal{A}_{2,1})x_{n+1/2} + (\mathcal{B}_{2,2} + \frac{1}{2}\mathcal{A}_{2,2})q_{n+1/2} \nonumber \\
&\quad - \int_{t_n}^{t_{n+1}} \theta(t/\varepsilon)(t-t_{n+1/2})dtJb_n, \nonumber \\
\log(r_{n+1}) - \log(r_n) &= \int_{t_n}^{t_{n+1}} b(s)ds \cdot (x_{n+1} - x_n), \tag{58}
\end{align}
$$

where $b_n = \nabla\phi(x_n)$ and $t_{n+1/2} = (t_n + t_{n+1})/2$. Let us remark that the only requirement for the approximation of $\int_{t_n}^{t_{n+1}} b(s)ds$ is to be second order accurate approximation of $\int_{t_n}^{t_{n+1}} \nabla\phi(x(s))ds$. Below, we present two ways to approximate $\int_{t_n}^{t_{n+1}} b(s)ds$.

**Choice 1** The first choice proposed by [31] is

$$
\frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} b(s) ds \approx \nabla\phi(x_{n+1/2}) = -\frac{1}{2}\nabla\phi(x_{n-1}) + \frac{3}{2}\nabla\phi(x_n). \quad (59)
$$

**Choice 2** The second choice follows the strategy we used to construct uniformly accurate scheme in the nonlinear case. First, we have

$$
\frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \nabla\phi(x(s))ds = \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} \nabla\phi(x(t_n) + h(s))ds,
$$

where $h(s) = \int_{t_n}^{s} q(s_1)ds_1 + \int_{t_n}^{s} \theta(s_1/\varepsilon)Jx(s_1)ds_1$ is obtained by integrating the equation on $x$ in (53). In view of a second order numerical scheme, we approximate $h(s)$ by $\tilde{h}(s)$ given by

$$
\tilde{h}(s) = (s - t_n)q(t_n) + \int_{t_n}^{s} \theta(s_1/\varepsilon)Jds_1x(t_n). \quad (60)
$$

Since $\phi$ is assumed to be smooth, we perform a Taylor expansion

$$
\begin{align*}
\int_{t_n}^{t_{n+1}} \nabla\phi(x(t_n) + h(s))ds &\approx \int_{t_n}^{t_{n+1}} \nabla\phi(x(t_n) + \tilde{h}(s))ds \\
&\approx \Delta t \nabla\phi(x(t_n)) + \nabla^2\phi(x(t_n)) \int_{t_n}^{t_{n+1}} \tilde{h}(s)ds \\
&\approx \Delta t \nabla\phi(x(t_n)) + \nabla^2\phi(x(t_n)) \int_{t_n}^{t_{n+1}} \left[ (s-t_n)q(t_n) + \int_{t_n}^{s} \theta(s_1/\varepsilon)Jds_1x(t_n) \right] ds \\
&\approx \Delta t \nabla\phi(x(t_n)) + \nabla^2\phi(x(t_n)) \left[ \frac{\Delta t^2}{2}q(t_n) + \int_{t_n}^{t_{n+1}} \int_{t_n}^{s} \theta(s_1/\varepsilon)Jds_1dsx(t_n) \right].
\end{align*}
$$
As a conclusion, $\frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} b(s)ds$ is chosen as follows

$$ \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} b(s)ds = \nabla\phi(x_n) + \frac{\Delta t}{2}\nabla^2\phi(x_n)q_n + \frac{1}{\Delta t}\nabla^2\phi(x_n)\int_{t_n}^{t_{n+1}}\int_{t_n}^{s}\theta(s_1/\epsilon)Jds_1dsx_n. \quad (61) $$

We end this part by identifying the limit of (58) when $\epsilon \to 0$ with the energy preserving scheme (48).

**Proposition 5.** The numerical scheme for charged particle model (58) with the approximation given by (61) degenerates when $\epsilon \to 0$ to the energy preserving numerical scheme for averaged model (48) with $\bar{b}$ given by (52).

*Proof.* The only point to check is the limit of (61) when $\epsilon \to 0$. From the definition of given by (61), we get (assuming $x_n \to \bar{x}_n$ and $q_n \to \bar{q}_n$ when $\epsilon \to 0$)

$$ \begin{aligned} \frac{1}{\Delta t} \int_{t_n}^{t_{n+1}} b(s)ds &= \nabla\phi(\bar{x}_n) + \nabla^2\phi(\bar{x}_n)\left[\frac{\Delta t}{2}\bar{q}_n + \frac{1}{\Delta t}\int_{t_n}^{t_{n+1}}(s-t_n)ds\langle\theta\rangle J\bar{x}_n\right] + \mathcal{O}(\epsilon) \\ &= \nabla\phi(\bar{x}_n) + \nabla^2\phi(\bar{x}_n)\frac{\Delta t}{2}\left[\bar{q}_n + \langle\theta\rangle J\bar{x}_n\right] + \mathcal{O}(\epsilon) = \bar{b} + \mathcal{O}(\epsilon), \end{aligned} $$

which is the definition (52) of $\bar{b}$ given for the SAV averaged model. $\square$

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

# 5 Particle-In-Cell

In this section, we are going to extend the numerical methods studied in previous sections for solving the Vlasov-Poisson equation (4) with a fast oscillating magnetic field. Indeed, using [12] and based on the two-scale convergence [1], the solution to (4) can be proved to converge to a Vlasov equation whose characteristics are given by (46). Hence, we extend our schemes developed in the ODE framework to the Vlasov case by considering the Particle-in-Cell (PIC) discretization. The PIC discretization approximates the unknown distribution $f^\epsilon(t, x, v)$ of (4) by a sum of Dirac masses located at $(x_k(t), v_k(t))$ with weight $\omega_k > 0$ for $k = 1, \dots, N_p$ and $N_p \in \mathbb{N}$ as

$$ f_p^\epsilon(t, x, v) = \sum_{k=1}^{N_p} \omega_k \delta(x - x_k(t)) \delta(v - v_k(t)), \quad t \ge 0, x, v \in \mathbb{R}^2. \quad (62) $$

The weights $\omega_k$ and initial values of the particles $x_{k,0}, v_{k,0}$ for $k = 0, \dots, N_p$ are prescribed according to the initial condition $f^\epsilon(t = 0, x, v) = f^{\text{in}}(x, v)$ in (4). Considering uniform weights, we have

$$ \omega_k = \frac{1}{N_p} \int_{\mathbb{R}^2 \times \mathbb{R}^2} f^{\text{in}}(x, v)dxdv, \quad k = 1, \dots, N_p. $$

Standard sampling techniques like the Monte Carlo type rejection sampling method found in standard textbooks [32] can be employed to obtain the initial data for position and velocity of the particles $x_{k,0}, v_{k,0}$ for $k = 1, \dots, N_p$. However, when $f^{\text{in}}$ takes the form $f^{\text{in}}(x, v) = \chi(x)M(v)$ with $\chi$ periodic and $M$ Maxwellian (which is a widely used case in plasma physics), specific deterministic techniques can be employed using the inversion of cumulative distribution function (see [10]).

The dynamics of macro-particles is then given by the characteristic equations (63) which are coupled to the Poisson equation

$$ -\Delta_x \phi^\epsilon(t, x) = \nabla_x \cdot E^\epsilon(t, x) = \rho^\epsilon(t, x) - 1, \quad \text{with } \rho^\epsilon(t, x) = \int_{\mathbb{R}^2} f^\epsilon(t, x, v)dv, $$

through the potential of electric field as $E^\epsilon = -\nabla_x \phi^\epsilon$. From the positions $\{x_k(t)\}_{k=1,\dots,N_p}$ of the particles at time $t > 0$, the density $\rho^\epsilon$ is approximated by $\rho_p$ given by

$$ \rho^\epsilon(t, x) \approx \rho_p^\epsilon(t, x) = \sum_{k=1}^{N_p} \omega_k \delta(x - x_k(t)) \quad x \in \mathbb{R}^2. $$

This enables to evaluate the density on a uniform grid so that the Poisson equation can be also solved on a mesh grid of $x$ in $\mathbb{R}^2$ to get $E_p^\epsilon(t, x) \approx E^\epsilon(t, x)$. Finally, an interpolation of $E_p^\epsilon(t, x)$ is required at each
particle position in the pusher. In practice, the Dirac function $\delta(x)$ is approximated by B-spline function $S^m(x)$ ($m \in \mathbb{N}$) given by [32]:

$$S^0(x) := \begin{cases} \frac{1}{\Delta x}, & |x| \le \frac{\Delta x}{2}, \\ 0, & \text{else,} \end{cases} \qquad S^m(x) := \frac{1}{\Delta x} \int_{x-\frac{\Delta x}{2}}^{x+\frac{\Delta x}{2}} S^{m-1}(y)dy, \quad m \ge 1.$$

The case in two dimensions is obtained by tensor product.

Under the PIC discretization (62), the characteristic equations of (4) read for $k = 1, \dots, N_p$

$$\dot{x}_k(t) = v_k(t), \qquad (63a)$$

$$\dot{v}_k(t) = E(t, x_k(t)) + \frac{B}{2\varepsilon}\theta'(t/\varepsilon)Jx_k(t) + B\theta(t/\varepsilon)Jv_k(t), \quad t > 0, \qquad (63b)$$

$$x_k(0) = x_{k,0}, \quad v_k(0) = v_{k,0}, \qquad (63c)$$

where $B$ denotes the magnetic field amplitude (the magnetic field $\mathbf{B} = \theta(t/\varepsilon)(0,0,B)$ is considered). Under suitable condition on $f^{\text{in}}$ (see [12]), we will consider the Vlasov equation on $g^{\varepsilon}(t,x,q) = f^{\varepsilon}(t,x,v)$ using the change of variables $q = v - B/2\theta(t/\varepsilon)Jx$. Hence, we consider the following ODE system (for $k = 1, \dots, N_p$)

$$\dot{x}_k(t) = q_k(t) + \frac{B}{2}\theta(t/\varepsilon)Jx_k(t), \qquad (64a)$$

$$\dot{q}_k(t) = E(t, x_k(t)) + \frac{B}{2}\theta(t/\varepsilon)Jq_k(t) + \frac{B^2}{4}\theta^2(t/\varepsilon)J^2x_k(t), \quad t > 0, \qquad (64b)$$

$$x_k(0) = x_{k,0}, \quad q_k(0) = v_{k,0} - \frac{B}{2}\theta(t/\varepsilon)Jx_{k,0}, \qquad (64c)$$

Hence, the PIC method is applied to $g^{\varepsilon}(t,x,q) \approx g_p(t,x,q) = \sum_{k=1}^{N_p} \omega_k \delta(x - x_k(t))\delta(q - q_k(t))$ where $x_k, q_k$ are solution to (64) for which we developed uniformly accurate numerical schemes in the previous section.

# 6 Numerical results

In this section, we present some numerical results to illustrate the explicit and midpoint uniformly accurate methods for the systems (41) and (8). We also present the numerical results for Vlasov-Poisson equation by employing Particle-In-Cell method along with our uniformly accurate methods for ODE systems.

## 6.1 ODE systems

The properties of numerical schemes for the ODE systems are validated in this subsection.

### UA property

The uniform accuracy (UA) properties of both explicit and energy preserving numerical schemes are presented in the linear case $\dot{U}(t) = A(t/\varepsilon)U(t)$, $U(t=0) = U_0, t \in [0, T], T = 1$. In Figures 1 to 3, we plot the $L^2$ error to illustrate the order of (uniform) accuracy. More precisely, in Figure 1, the first order explicit UA scheme is considered with $U(t) = (x(t), q(t)) \in \mathbb{R}^4$, $\theta(s) = 1 + \cos(s)$ in the matrix $A \in \mathcal{M}_{4,4}(\mathbb{R})$ given by

$$A(t/\varepsilon) = \begin{bmatrix} \theta(t/\varepsilon)J & 1 \\ \theta^2(t/\varepsilon)J^2 & \theta(t/\varepsilon)J \end{bmatrix}. \qquad (65)$$

In Figure 2, the standard midpoint and UA midpoint schemes are considered with the same ODE as before except $\theta(s) = \cos(s)$. It is observed that second order accuracy is obtained for both schemes. Finally, fourth order UA scheme is considered in Figure 3 for the scalar model $\dot{U}(t) = [2 + 0.5 \cos^2(t/\varepsilon)]U(t)$. For all the UA schemes studied in Figures 1 to 3, the uniform accuracy is recovered.

In Figures 4,5,6, the nonlinear case is considered: $\dot{U}(t) = A(t/\varepsilon)U(t) + g(U(t))$, $U(t=0) = U_0, t \in [0, T], T = 1$, with $U(t) = (x(t), q(t)) \in \mathbb{R}^4$ and $g(U(t)) = (0, 0, \cos x_1 \sin x_2 + x_1 + x_1^3, \sin x_1 \cos x_2 + x_2 + x_2^3)^T \in \mathbb{R}^4$. Finally, the matrix $A$ given by (65) is considered with $\theta(s) = \cos(s)$. Figure 4 depicts the results obtained by the second order explicit UA scheme (43). The expected second order uniform accuracy is recovered. Figures 5,6 correspond to midpoint second order scheme for the SAV reformulated system given by (58). We compare the two choices of $b$: in Figure 5, the results are obtained with the first choice for the approximation of $\int_{t_n}^{t_{n+1}} b(s)ds$ given by (59) whereas Figure 6 involves the approximation given by (61). It is observed that both choices reach the second order accuracy. Further, we also observe that the numerical scheme (48) for the asymptotic model with both choices of $b$ (given by (50) or (51) preserve energy up to machine accuracy for long time (figures not shown).
## Geometric property of the averaged model

In this part, we present the phase plots and frequency spectra obtained by energy preserving methods for the approximation of the following models: $\dot{U} = A(t/\epsilon)U$, A is given by (65) with $\theta$ replaced by $B\theta$, $B \in \mathbb{R}$ being the magnitude of magnetic field. When $\epsilon$ goes to zero, the asymptotic model is determined by the eigenvalues of $\langle A \rangle$ which can be computed for a given function $\theta$.

Figure 7 corresponds to an explicit first order uniformly accurate scheme for the ODE system with $B = 1$, $\Delta t = 0.001$, $\theta(s) = 1 + \cos s$, and the final time is $T = 100$. We consider $\epsilon = 0.1$ and display the time history of $x_1$ together with its averaged counterpart $\bar{x}_1$ (see Figure 7-(a)), the phase plot of $(x_1, x_2)$ (see Figure 7-(b)) and the phase plot of $(\bar{x}_1, \bar{x}_2)$ (see Figure 7-(c)). In Figure 7-(d) (resp. (e)), we display the discrete Fourier transform of the sequence $(x_1)_n$ (resp. $(\bar{x}_1)_n$). Two main peaks can be observed corresponding to the main periods of the numerical solutions which are in very good agreement with the frequencies computed from the eigenvalues of $\langle A \rangle$ in the averaged model (equal to $i0.224$ and $i2.224$, corresponding to the arrows). Let us remark that additional frequencies are observed for the spectrum of $(x_1)_n$ due to oscillations in $\epsilon$. Figures 8 to 10 show numerical results for the nonlinear problem obtained with the second order uniformly accurate energy preserving scheme with $\Delta t = 0.1$, $\theta(s) = \cos s$, for different values of $B = 0.5, 1$ and $5$ respectively. Different values of $\epsilon = 0.1, 0.001$ are considered to study the behaviour of the system. Note that for this choice of $\theta$, there is only one frequency corresponding to the same imaginary value present in all the eigenvalues of $\langle A \rangle$ (equal to $i0.707107B$ which is well captured by the numerical scheme, not shown). In Figures 8 to 10, the phase plots of $(x_1, x_2)$ and $(\bar{x}_1, \bar{x}_2)$ are displayed to investigate the confinement property with respect to both $\epsilon$ and $B$. We observe the domain of $(x_1, x_2)$ shrinks with increase in $B$ which illustrates an improved confinement. When $\epsilon$ decreases, the perturbations in $(x_1, x_2)$ become smaller, resulting in better confinement. Further, energy is preserved up to machine accuracy (figures not shown) for the scheme (48) of the averaged model.

## 6.2 Vlasov-Poisson solver: Particle-In-Cell method

In this section, we present some numerical results obtained for the Vlasov-Poisson equation by using PIC method that employs our uniformly accurate numerical methods for solving the ODE system (64). We consider the following initial condition

$$f^{\text{in}}(x, v) = \left(1 + \xi_1 \cos(k_1 x_1)\right) \left(1 + \xi_2 \cos(k_2 x_2)\right) \frac{1}{2\pi} e^{-|v|^2/2}. \quad (66)$$

We also note that when $B = 0$ in (64), the problem degenerates to the Vlasov-Poisson equation. Thus, in this case, the initial condition (66) corresponds to a Landau damping test which enables to validate our code by comparing our results with those in literature. Then, we investigate how the Landau damping behaves for non-zero values of $B$. We thus plot the time history electric energy $\int \int |E|^2 dx_1 dx_2$ (semi-log scale). In Figures 11 and 12, we consider the one dimensional case with $\xi_2 = 0$ and $k_2 = 2\pi$, whereas $\xi_1 = 0.05$. The domain is $x_1 \in [0, 2\pi/k_1]$, $x_2 \in [0, 2\pi/k_2]$. The number of points along $x_1$ (resp. $x_2$) direction is $N_1 = 128$ (resp. $N_2 = 4$), the number of particles is $N_p = 100N_1N_2$, and time step is $\Delta t = 0.01$. Figure 11 depicts the time history electric energy with $B = 0$ for different values of $k_1$ ($k_1 = 0.5, 0.4, 0.3$). The electric energy is expected to decrease exponentially fast in time with a rate which can be computed from the linear theory (see [32]). In Figure 11, it is observed that the obtained slopes are close to the theoretical values of damping observed in [32]. In Figure 12, we consider $k_1 = 0.3$ case for non-zero values of $B = 0.01, 0.05, 0.1$ and $0.15$, with magnetic field's oscillatory parameter $\epsilon$ taken as 0.001. When $B$ becomes larger $B \ge 0.1$, the damping disintegrates.

In Figures 13 and 14, we consider the two dimensional case with $\xi_1 = \xi_2 = 0.05$. The domain is $(x_1, x_2) \in [0, 2\pi/k_1] \times [0, 2\pi/k_2]$ (the numerical parameters are $N_1 = N_2 = 128$, $N_p = 100N_1N_2$ and $\Delta t = 0.01$). Figure 13 depicts the Landau damping phenomena with $B = 0$ for different values of $k_1 = k_2$ such as $0.5, 0.4, 0.3$. There is not too much references for two-dimensional Landau damping test but we observe that the rate obtained with $k_1 = k_2 = 0.5$ is very close to the theoretical one $-0.15$ (see [27]).

In fig. 14, we consider the $k_1 = k_2 = 0.3$ case for non-zero values of $B = 0.01, 0.05, 0.1$ and $0.15$, with magnetic field's oscillatory parameter $\epsilon$ taken as 0.001. Similar to the one dimensional case, the electric energy enjoys a very different behavior when $B$ increases.

## 7 Conclusion

In this work, we design and analyse uniformly accurate numerical schemes for highly oscillatory differential equations of the form $\dot{U} = A(t/\epsilon)U + g(U)$. Both explicit and linearly implicit schemes are studied, according
<figure><img src="image_3.png" /><figcaption>(a) $L^2$ error of $x_1, x_2$ as a function of $\Delta t$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) $L^2$ error of $q_1, q_2$ as a function of $\Delta t$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(d) $L^2$ error of $q_1, q_2$ as a function of $\varepsilon$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(a) error with $\Delta t$ - midpoint</figcaption></figure>

<caption>(b) error with $\Delta t$ - UA midpoint</caption>

<caption>(c) error with $\varepsilon$ - UA midpoint</caption>

<figure><img src="image_7.png" /><figcaption>(a) $L^2$ error of $y$ as a function of $\Delta t$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) $L^2$ error of $y$ as a function of $\varepsilon$</figcaption></figure>

<caption>Figure 3: Linear scalar case, fourth order explicit scheme: uniform accuracy with $\Delta t$.</caption>
<figure><img src="image_9.png" /><figcaption>Figure 4: Nonlinear case, second order explicit scheme: Second order uniform accuracy with $\Delta t$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 5: Nonlinear case, midpoint SAV-scheme (first choice for $b$): Second order uniform accuracy with $\Delta t$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 6: Nonlinear case, midpoint SAV-scheme (second choice for $b$): Second order uniform accuracy with $\Delta t$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 7: Explicit first order uniformly accurate scheme - There are two different frequencies in the sequence $\bar{X}_1$ (asymptotic model); the $\varepsilon$-model also has small frequencies due to perturbations arising from $\varepsilon$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 8: $B = 0.5$, phase plots of $(x_1, x_2)$ for $\varepsilon = 0.1, 0.001$ and of $(\bar{X}_1, \bar{X}_2)$ obtained with the energy preserving second order uniformly accurate scheme for $\dot{U} = A(t/\varepsilon)U$ with $A$ given by (65) and $\theta(s) = \cos(s)$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9: $B = 1$, phase plots of $(x_1, x_2)$ for $\varepsilon = 0.1, 0.001$ and of $(\bar{X}_1, \bar{X}_2)$ obtained with the energy preserving second order uniformly accurate scheme for $\dot{U} = A(t/\varepsilon)U$ with $A$ given by (65) and $\theta(s) = \cos(s)$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 10: $B = 5$, phase plots of $(x_1, x_2)$ for $\epsilon = 0.1, 0.001$ and of $(\bar{X}_1, \bar{X}_2)$ obtained with the energy preserving second order uniformly accurate scheme for $\dot{U} = A(t/\epsilon)U$ with $A$ given by (65) and $\theta(s) = \cos(s)$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 11: Vlasov-Poisson (4) with $B = 0$. Landau damping test ($\xi_1 = 0.05, \xi_2 = 0$): time history of the electric energy for $k_1 = 0.5$ (left), $k_1 = 0.4$ (middle) and $k_1 = 0.3$ (right).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 12: Vlasov-Poisson (4) with $\epsilon = 0.001$. Landau damping test ($\xi_1 = 0.05, \xi_2 = 0, k_1 = 0.3$): time history of the electric energy for different values of $B$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 13: Vlasov-Poisson (4) with $B = 0$. Landau damping test ($\xi_1 = \xi_2 = 0.05$): time history of the electric energy for $k_1 = k_2 = 0.5$ (left), $k_1 = k_2 = 0.4$ (middle) and $k_1 = k_2 = 0.3$ (right).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 14: Vlasov-Poisson (4) with $\epsilon = 0.001$. Landau damping test ($\xi_1 = \xi_2 = 0.05$, $k_1 = k_2 = 0.3$): time history of the electric energy for different values of $B$.</figcaption></figure>

to the structure of the equation. In particular, when the averaged model enjoys a Hamiltonian structure, midpoint type schemes are proposed which degenerate into energy preserving schemes for the averaged model. Several numerical results illustrate the accuracy and geometric properties of the scheme. Further, the solvers have been used within PIC framework to solve the Vlasov-Poisson equation.

Several extensions to this work can be performed: for instance, the application to the dynamics of charged particles under the influence of strong magnetic fields, but also, links with Magnus integrators will be investigated.

## A First order uniform accuracy of the standard midpoint scheme

In this first Appendix, we prove that the standard midpoint scheme (20) is only first order UA.

**Theorem 6.** The numerical scheme (20) is first order uniformly accurate.

*Proof.* Considering the error equation by subtracting (23) with (20), we obtain

$$ 
\begin{aligned}
e_{n+1} - e_n &= \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left(U(s) - \frac{U_n + U_{n+1}}{2}\right) ds \\
&= \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \left(\frac{e_n + e_{n+1}}{2}\right) + \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left(\frac{U(s) - U(t_n)}{2} + \frac{U(s) - U(t_{n+1})}{2}\right) ds \\
&= \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \left(\frac{e_n + e_{n+1}}{2}\right) + \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left(\int_{t_n}^{s} A(s_1/\epsilon) U(s_1) ds_1 - \int_{s}^{t_{n+1}} A(s_1/\epsilon) U(s_1) ds_1\right) ds. \quad (67)
\end{aligned}
$$

Then, as the solution $U$ is uniformly bounded according to lemma B-1, the norm of $e_{n+1}$ becomes,

$$ \|e_{n+1}\| \le \frac{1 + C \Delta t}{1 - C \Delta t} \|e_n\| + \frac{C}{1 - C \Delta t} \Delta t^2. $$

We then conclude by Gronwall lemma that $\|e_n\| \le C \Delta t$, with $C$ independent of $\Delta t$ and $\epsilon$.

However, it is well known that the midpoint scheme is second order for non stiff problems ($A$ independent of time or for $\epsilon$ large). Let us thus focus on the last term in (67). It is possible to replace $U(s_1)$ by $U(t_n) + \int_{t_n}^{s_1} A(s_2/\epsilon)U(s_2)ds_2$ to get

$$ 
\begin{aligned}
e_{n+1} - e_n &= \int_{t_n}^{t_{n+1}} A(s/\epsilon) ds \left(\frac{e_n + e_{n+1}}{2}\right) + \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left(\int_{t_n}^{s} A(s_1/\epsilon) ds_1 - \int_{s}^{t_{n+1}} A(s_1/\epsilon) ds_1\right) ds U(t_n) \\
&\quad + \frac{1}{2} \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left(\int_{t_n}^{s} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2 ds_1 - \int_{s}^{t_{n+1}} A(s_1/\epsilon) \int_{t_n}^{s_1} A(s_2/\epsilon) U(s_2) ds_2 ds_1\right) ds.
\end{aligned}
$$

The last term is of order $\Delta t^3$ uniformly in $\epsilon$ due the uniform boundedness of $U$ and of $A$. The loss of uniform accuracy thus should come from the second term. A direct analysis gives a bound $\Delta t^2$ uniformly in $\epsilon$, but this does not take the non stiff case into account. A Taylor expansion of $A(s_1/\epsilon)$ gives

$$ A(s_1/\epsilon) = A(t_{n+1/2}/\epsilon) + \frac{(s_1 - t_{n+1/2})}{\epsilon} A'(\xi/\epsilon) $$

so that

$$ 
\begin{aligned}
\int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{t_n}^{s} A(s_1/\epsilon) ds_1 ds &= \int_{t_n}^{t_{n+1}} A(s/\epsilon)(s - t_n) A(t_{n+1/2}/\epsilon) ds + \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{t_n}^{s} \frac{(s_1 - t_{n+1/2})}{\epsilon} A'(\xi/\epsilon) ds_1 ds \\
\int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{s}^{t_{n+1}} A(s_1/\epsilon) ds_1 ds &= \int_{t_n}^{t_{n+1}} A(s/\epsilon)(t_{n+1} - s) A(t_{n+1/2}/\epsilon) ds + \int_{t_n}^{t_{n+1}} A(s/\epsilon) \int_{s}^{t_{n+1}} \frac{(s_1 - t_{n+1/2})}{\epsilon} A'(\xi/\epsilon) ds_1 ds
\end{aligned}
$$
and the difference becomes

$$ \left| \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left[ \int_{t_n}^{s} A(s_1/\epsilon) ds_1 - \int_{s}^{t_{n+1}} A(s_1/\epsilon) ds_1 \right] ds \right| \le 2 \left| \int_{t_n}^{t_{n+1}} A(s/\epsilon)(s - t_{n+1/2}) A(t_{n+1/2}/\epsilon) ds \right| + \left| \int_{t_n}^{t_{n+1}} A(s/\epsilon) \left[ \int_{t_n}^{s} \frac{(s_1 - t_{n+1/2})}{\epsilon} A'(\xi/\epsilon) ds_1 - \int_{s}^{t_{n+1}} \frac{(s_1 - t_{n+1/2})}{\epsilon} A'(\xi/\epsilon) ds_1 \right] ds \right|. $$

The last term is bounded by $C\Delta t^3/\epsilon$. Let us now focus on the first term on the right hand side:

$$ \begin{aligned} \left| \int_{t_n}^{t_{n+1}} A(s/\epsilon)(s - t_{n+1/2})ds \right| &= \left| \int_{t_n}^{t_{n+1}} \left[ A(t_{n+1/2}/\epsilon) + \frac{(s - t_{n+1/2})}{\epsilon} A'(\xi/\epsilon) \right] (s - t_{n+1/2})ds \right| \\ &= \left| \int_{t_n}^{t_{n+1}} \frac{(s - t_{n+1/2})^2}{\epsilon} A'(\xi/\epsilon)ds \right| \le C \frac{\Delta t^3}{\epsilon}. \end{aligned} $$

We thus conclude with the following relation on the error

$$ \|e_{n+1}\| \le \frac{1 + C\Delta t}{1 - C\Delta t} \|e_n\| + \frac{C\Delta t^3/\epsilon}{1 - C\Delta t}, $$

from which we get, using Gronwall lemma $\|e_n\| \le C\Delta t^2/\epsilon$.

Gathering the two estimates $\|e_n\| \le C\Delta t$ and $\|e_n\| \le C\Delta t^2/\epsilon$ leads to $\|e_n\| \le C \min(\Delta t, \Delta t^2/\epsilon)$ and we conclude that the midpoint method is first order uniformly accurate. $\square$

## B Technical lemmas

In this Appendix, we gather several technical lemmas which are used several times in different proofs. First, a very standard lemma ensures that the solution of the ODE is uniformly bounded.

**Lemma 1** (Uniform boundedness (with $\epsilon$) of the exact solution). Consider $U(t)$ the exact solution to (11), where $g(t, U(t))$ is a Lipschitz continuous function with Lipschitz constant $K$ and $A$ is a bounded matrix-valued function. Then there exist positive constants $K_1$ and $K_2$, independent of $\epsilon$, such that for $t \ge 0$, $\|U(t)\| \le K_1 \exp \left( \int_0^t K_2 ds \right)$.

The proof is classical and derives from the Gronwall lemma and the following representation formula:

$$ U(t) = U(0) + \int_{0}^{t} A(s/\epsilon) U(s) ds + \int_{0}^{t} g(s, U(s)) ds. $$

The asymptotic expansions of averaged value for highly oscillating functions is useful.

**Lemma 2.** Let $\theta$ be a $P$-periodic smooth function and set $\langle\theta\rangle = \frac{1}{P} \int_0^P \theta(s) ds$ as its average. Then the following asymptotic expansions are satisfied ($t_n = n\Delta t$ for any integer $n$):

$$ \int_a^b \theta(t/\epsilon) dt = (b-a)\langle\theta\rangle + \mathcal{O}((b-a)\epsilon), \quad (68) $$

$$ \int_{t_n}^{t_{n+1}} \int_{t_n}^{t} \theta(t_1/\epsilon) dt_1 dt = \frac{\Delta t^2}{2} \langle\theta\rangle + \mathcal{O}(\Delta t^2 \epsilon), \quad (69) $$

$$ \int_{t_n}^{t_{n+1}} \int_{t}^{t_{n+1}} \theta(t_1/\epsilon) dt_1 dt = \frac{\Delta t^2}{2} \langle\theta\rangle + \mathcal{O}(\Delta t^2 \epsilon), \quad (70) $$

$$ \int_{t_n}^{t_{n+1}} \theta(t/\epsilon)(t - t_n) dt = \frac{\Delta t^2}{2} \langle\theta\rangle + \mathcal{O}(\Delta t \epsilon), \quad (71) $$

$$ \int_{t_n}^{t_{n+1}} \theta(t/\epsilon)(t_{n+1} - t) dt = \frac{\Delta t^2}{2} \langle\theta\rangle + \mathcal{O}(\Delta t \epsilon). \quad (72) $$

*Proof.* Without loss of generality, we may assume the period $P = 2\pi$ and consider then the Fourier series expansion of $\theta(s)$, with coefficients $C_k = \frac{1}{2\pi} \int_0^{2\pi} \theta(s)e^{-iks} ds$, $k \in \mathbb{Z}$, having appropriate summability properties inherited from the smoothness of $\theta$:

$$ \theta(s) = \sum_{k \in \mathbb{Z}} C_k e^{iks} = \langle\theta\rangle + \sum_{k \in \mathbb{Z} \setminus \{0\}} C_k e^{iks}, $$
we obtain

$$ \int_a^b (\theta(t/\epsilon) - \langle\theta\rangle) dt = \sum_{k \in \mathbb{Z} \setminus \{0\}} C_k \int_a^b e^{ikt/\epsilon} dt = \epsilon \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{1}{ik} C_k [e^{ikb/\epsilon} - e^{ika/\epsilon}] = \mathcal{O}((b-a)\epsilon). $$

Thus, the identities (68), (69) and (70) become evident:

$$ \begin{aligned} \int_a^b \theta(t/\epsilon) dt &= (b-a)\langle\theta\rangle + \mathcal{O}((b-a)\epsilon), \\ \int_{t_n}^{t_{n+1}} \int_{t_n}^t \theta(t_1/\epsilon) dt_1 dt &= \int_{t_n}^{t_{n+1}} [\langle\theta\rangle(t-t_n) + \mathcal{O}(\Delta t \epsilon)] dt = \frac{\Delta t^2}{2} \langle\theta\rangle + \mathcal{O}(\Delta t^2 \epsilon), \\ \int_{t_n}^{t_{n+1}} \int_{t_n}^{t_{n+1}} \theta(t_1/\epsilon) dt_1 dt &= \int_{t_n}^{t_{n+1}} [\langle\theta\rangle(t_{n+1}-t) + \mathcal{O}(\Delta t \epsilon)] dt = \frac{\Delta t^2}{2} \langle\theta\rangle + \mathcal{O}(\Delta t^2 \epsilon). \end{aligned} $$

To prove (71) and (72), we use similar calculations in the following way

$$ \begin{aligned} \int_{t_n}^{t_{n+1}} (\theta(t/\epsilon) - \langle\theta\rangle)(t-t_n)dt &= \int_{t_n}^{t_{n+1}} \sum_{k \in \mathbb{Z} \setminus \{0\}} C_k e^{ikt/\epsilon} (t-t_n)dt \\ &= \epsilon \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{1}{ik} C_k [e^{ikt/\epsilon}(t-t_n)]_{t_n}^{t_{n+1}} - \epsilon \int_{t_n}^{t_{n+1}} \sum_{k \in \mathbb{Z} \setminus \{0\}} \frac{1}{ik} C_k e^{ikt/\epsilon} dt = \mathcal{O}(\epsilon \Delta t). \end{aligned} $$

□

## References

[1] ALLAIRE, G. Homogenization and two-scale convergence. SIAM J. Math. Anal. 23 (1992), 1482–1518.

[2] ASTUTO, C. High order multiscale methods for advection-diffusion equation with highly oscillatory boundary condition, 2024.

[3] ASTUTO, C., LEMOU, M., AND RUSSO, G. Time multiscale modeling of sorption kinetics i: uniformly accurate schemes for highly oscillatory advection-diffusion, 2023.

[4] BAO, W., AND CAI, Y. Uniform and optimal error estimates of an exponential wave integrator sine pseudospectral method for the nonlinear Schrödinger equation with wave operator. SIAM J. Numer. Anal. 52 (2014), 1103–1127.

[5] BAO, W., CAI, Y., JIA, X., AND TANG, Q. A uniformly accurate multiscale time integrator pseudospectral method for the Dirac equation in the nonrelativistic limit regime. SIAM J. Numer. Anal. 54 (2016), 1785–1812.

[6] BAO, W., CAI, Y., AND ZHAO, X. A uniformly accurate multiscale time integrator pseudospectral method for the Klein–Gordon equation in the nonrelativistic limit regime. SIAM J. Numer. Anal. 52 (2014), 2488–2511.

[7] BAO, W., DONG, X., AND ZHAO, X. An exponential wave integrator sine pseudospectral method for the Klein–Gordon–Zakharov system. SIAM J. Scientific Comput. 35 (2013), 2903–2927.

[8] BAO, W., DONG, X., AND ZHAO, X. Uniformly correct multiscale time integrators for highly oscillatory second order differential equations. J. Math. Study 47 (2014).

[9] BAUMSTARK, S., FAOU, E., AND SCHRATZ, K. Uniformly accurate exponential-type integrators for Klein-Gordon equations with asymptotic convergence to the classical NLS splitting. Math. Comp. 87 (2018).

[10] BIRDSALL, C., AND LANGDON, A. Plasma Physics via Computer Simulation. Adam Hilger, 1991.

[11] BOGOLIUBOV, N. N., AND MITROPOLSKY, Y. A. Asymptotic methods in the theory of non-linear oscillations. Gordon and Breach Science Publishers, Inc., New York, 1961.
[12] BOSTAN, M. Transport of Charged Particles Under Fast Oscillating Magnetic Fields. *SIAM J. Math. Anal.* **44**, 3 (2012), 1415–1447.

[13] CAI, Y., AND GUO, Y. Uniformly accurate nested Picard integrators for a system of oscillatory ordinary differential equations. *BIT Numer Math* **61** (2021).

[14] CAI, Y., AND WANG, Y. Uniformly accurate nested Picard iterative integrators for the Dirac equation in the nonrelativistic limit regime. *SIAM J. Numer. Anal.* **57** (2019).

[15] CHARTIER, P., CROUSEILLES, N., LEMOU, M., AND MÉHATS, F. Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations. *Numer. Math.* **129** (2015).

[16] CHARTIER, P., CROUSEILLES, N., LEMOU, M., MÉHATS, F., AND ZHAO, X. Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field. *Math. Comp.* **88** (2019).

[17] CHARTIER, P., CROUSEILLES, N., LEMOU, M., MÉHATS, F., AND ZHAO, X. Uniformly accurate methods for three dimensional Vlasov equations under strong magnetic field with varying direction. *SIAM J. Sc. Comput.* **42** (2020).

[18] CHARTIER, P., CROUSEILLES, N., AND ZHAO, X. Numerical methods for the two-dimensional Vlasov-Poisson models in the finite Larmor radius approximation regime. *J. Comput. Phys.* **375** (2018).

[19] CHARTIER, P., LEMOU, M., MÉHATS, F., AND VILMART, G. A new class of uniformly accurate numerical schemes for highly oscillatory evolution equations. *FoCM* **20** (2020).

[20] CHARTIER, P., LEMOU, M., MÉHATS, F., AND ZHAO, X. Derivative-free high-order uniformly accurate schemes for highly oscillatory systems. *IMA J. Numer. Anal.* **42** (2022).

[21] CRESTETTO, A., CROUSEILLES, N., AND PREL, D. Multiscale numerical schemes for the collisional Vlasov equation in the finite Larmor radius approximation regime. *Multiscale Modeling and Simulation* **21**(3) (2023).

[22] CROUSEILLES, N., HIRSTOAGA, S., AND ZHAO, X. Multiscale Particle-In-Cell methods and comparisons for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field. *Comput. Phys. Comm.* **222** (2018).

[23] CROUSEILLES, N., JIN, S., LEMOU, M., AND MÉHATS, F. A micro-macro method for a kinetic graphene model in one-space dimension. *Multiscale Modeling and Simulation* **18** (2020).

[24] CROUSEILLES, N., LEMOU, M., AND MÉHATS, F. Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations. *J. Comput. Phys.* **248** (2013).

[25] CROUSEILLES, N., LEMOU, M., MÉHATS, F., AND ZHAO, X. Uniformly accurate forward semi-lagrangian methods for highly oscillatory Vlasov-Poisson equations. *Multiscale Modeling and Simulation* **15** (2017).

[26] CROUSEILLES, N., LEMOU, M., MÉHATS, F., AND ZHAO, X. Uniformly accurate Particle-In-Cell method for the long time solution of the two-dimensional Vlasov-Poisson equation with uniform strong magnetic field. *J. Comput. Phys.* **346** (2017).

[27] FILBET, F., AND SONNENDRUCKER, E. Comparison of eulerian vlasov solvers. *Comput. Phys. Comm.* **150** (2003), 247–266.

[28] HAIRER, E., LUBICH, C., AND WANNER, G. *Geometrical Numerical Analysis*. Springer, 2004.

[29] ISERLES, A., MUNTHE-KAAS, H., NORSETT, S., AND ZANNA, A. Lie-group methods. *Acta Numerica* **9** (2000), 215–365.

[30] MOCQUARD, Y., NAVARO, P., AND CROUSEILLES, N. *Hoodesolver.jl*: A julia package for highly oscillatory problems. *JOSS* **61**(6) (2021).

[31] SHEN, J., AND XU, J. Convergence and error analysis for the scalar auxiliary variable (sav) schemes to gradient flows. *SIAM J. Numer. Anal.* **56**, 5 (2018), 2895–2912.

[32] SONNENDRÜCKER, E. *Numerical Methods for Vlasov-Maxwell Equations*. Lecture notes, 2016.

[33] WANG, B., AND ZHAO, X. Geometric two-scale integrators for highly oscillatory system: uniform accuracy and near conservations. *SIAM J. Numer. Anal.* (2023).