# First-order two-scale simulations for charged particle dynamics in inhomogeneous strong magnetic fields

Emmanuel Frénod, Sever Adrian Hirstoaga

► To cite this version:

Emmanuel Frénod, Sever Adrian Hirstoaga. First-order two-scale simulations for charged particle dynamics in inhomogeneous strong magnetic fields. 2025. hal-05028088

HAL Id: hal-05028088

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

Preprint submitted on 9 Apr 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# First-order two-scale simulations for charged particle dynamics in inhomogeneous strong magnetic fields

Emmanuel Frénod *

Sever A. Hirstoaga †

April 9, 2025

## Abstract

We develop a multi-scale numerical approach to solve the Newton-Lorentz equations for the dynamics of a charged particle in inhomogeneous in space and time electromagnetic fields with a strong magnetic component. Based on the general setting of [11], we consider zero-order and first-order two-scale models derived from asymptotic expansions, to approximate the initial model. Then, we implement a two-scale numerical method as follows: first, with the help of symbolic software we obtain explicit forms of the right-hand side terms in the two-scale equations, and second, we classically use a Runge-Kutta integrator to solve the equations. The numerical experiments demonstrate the relevance of the first-order two-scale solution for various simulations of multi-scale plasma. Indeed, the first-order solution proves to be highly accurate when approximating the original model and computationally efficient thanks to its large time step compared to the cyclotron period.

**Keywords:** multi-scale model, Newton-Lorentz equation, first-order approximation, symbolic computation, Runge-Kutta integrator

## 1 Introduction

The motivation of this work is to efficiently numerically integrate a stiff Vlasov-Poisson system by means of a particle-in-cell method. While the kinetic equation models the evolution of a plasma under the influence of a strong magnetic field, the particle-in-cell algorithm represents the physical particles distribution by numerical particles which are tracked in time in the phase space [3, 17, 15]. More precisely, our general goal is to solve equations of the form

$$ \frac{\partial f_{\varepsilon}}{\partial t} + \mathbf{v} \cdot \nabla_x f_{\varepsilon} + \left( \mathbf{E}(t, \mathbf{x}) + \mathbf{v} \times \left( \mathbf{B}(\mathbf{x}) + \frac{1}{\varepsilon} m(t, \mathbf{x}) \right) \right) \cdot \nabla_{\mathbf{v}} f_{\varepsilon} = 0, \quad f_{\varepsilon}(t = s, \mathbf{x}, \mathbf{v}) = f_0(\mathbf{x}, \mathbf{v}), \quad (1) $$

where $0 < \varepsilon \ll 1$, $f_{\varepsilon} \equiv f_{\varepsilon}(t, \mathbf{x}, \mathbf{v})$ is the distribution function of a particle species (electrons or ions), $t$ stands for the time variable, $(\mathbf{x}, \mathbf{v})$ for the position and velocity variables, and $f_0$ is the initial condition. In addition, $(\mathbf{m}, \mathbf{E}) : \mathbb{R}^+ \times \mathbb{R}^3 \to \mathbb{R}^6$ and $\mathbf{B} : \mathbb{R}^3 \to \mathbb{R}^3$ are given vector fields

<footnote>* Université Bretagne Sud, Laboratoire de Mathématiques de Bretagne Atlantique, UMR 6205, F-56000 Vannes, France.</footnote>

<footnote>† Inria, Sorbonne Université, CNRS, Université de Paris, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France.</footnote>
such that $m$ is of constant norm 1. The parameter $1/\epsilon$ in front of the $m$ term refers to the fact that this magnetic term is assumed high with respect to the other terms. If we also consider the electric field generated by the particles, an additional term appears in the $E$ term above, calculated from the Poisson equation. The difficulty when solving (1) is that the solution contains both high-frequency and low-frequency time oscillations. Therefore, in order to cover both time scales, a particle-in-cell simulation should consist of a very large number of time iterations which results in a high computational load; see for instance [7, 14, 2, 9, 10] for different high performance computing or numerical strategies to deal with this issue. An alternative to avoid this situation is to develop methods or models that are free of the cyclotron frequency $1/\epsilon$. One strategy in this direction is to use asymptotic expansions in the small parameter $\epsilon$ to approximate the unknown $f_\epsilon$ by a sum in powers of $\epsilon$ with coefficients to be determined. For instance, gyrokinetic models can be built with this approach, by using the geometrical writing of the Vlasov equation and sophisticated tools as differential forms on symplectic manifolds (see for example [21, 19, 15, 5, 13]). We mention also the importance of determining orders 1 or 2 in these expansions to capture an oscillating part of $f_\epsilon$ which might play a role in its long-time behaviour (see [4, 6, 12, 8] among the large literature on the subject). A second strategy, which motivates this work, is as follows: when solving the Vlasov equation with the particle-in-cell method we are led to consider the characteristics of the equation. Then, we use asymptotic expansions for the unknowns of these ODEs in order to derive reduced models that are much less expensive than the initial full model. Finally, we can recover from these characteristics the solution $f_\epsilon$ to the Vlasov equation (see [11]).

In this framework, in the present paper we are solving numerically the following six dimensional dynamical system

$$ \left\{ \begin{aligned} \frac{d\mathbf{x}_\epsilon}{dt} &= \mathbf{v}_\epsilon, & \mathbf{x}_\epsilon(s) &= \mathbf{x}, \\ \frac{d\mathbf{v}_\epsilon}{dt} &= \mathbf{v}_\epsilon \times \left( \frac{1}{\epsilon} m(t, \mathbf{x}_\epsilon) + \mathbf{B}(\mathbf{x}_\epsilon) \right) + \mathbf{E}(t, \mathbf{x}_\epsilon), & \mathbf{v}_\epsilon(s) &= \mathbf{v}, \end{aligned} \right. \quad (2) $$

where the parameter $0 < \epsilon \ll 1$ and the vector fields $(m, E, B)$ are as above. The vector $(\mathbf{x}, \mathbf{v})$ is the initial condition given at the initial time $t = s \ge 0$. The system in (2) models the dynamics of a single charged particle in a strong magnetic field $\frac{1}{\epsilon}m$ and a weaker electromagnetic field $(E, B)$. The charge and the mass of the particle are both fixed to 1. In this context, $\mathbf{x}_\epsilon : \mathbb{R}^+ \to \mathbb{R}^3$ stands for the position unknown and $\mathbf{v}_\epsilon : \mathbb{R}^+ \to \mathbb{R}^3$ for the velocity unknown. As explained above, the multi-scale nature of the problem stems from the fact that the large magnetic field term introduces a rapid time scale, whereas the electric field term drives a much slower motion. Thus, any standard numerical integrator of (2) needs to use time steps that are of the order of the fastest oscillation if high accuracy is required, which leads to expensive simulations; the literature on this subject is extensive, see for example [18, 16, 20, 22, 23, 25]. Our goal is to propose a new numerical method whose time step is not constrained by $\epsilon$.

To this end, we rely the two-scale asymptotic expansions in the parameter $\epsilon$ of the unknowns in (2) to derive reduced models, whose importance lies in the fact that they are free from high frequency oscillations and approximate the model in (2) when $\epsilon$ vanishes. This approach was developed in [11], where abstract general formulae have been obtained for two-scale models to any order in $\epsilon$. In addition, theoretical results of convergence have been proven for these models to justify the approximations. In this context, the aim of the present work is as follows: first, to
identify illustrative test cases of strong magnetic fields relevant to physics, then to deduce the forms of the two-scale models of order 0 and 1 in these cases, and finally to solve them numerically by a suitable scheme in order to underline their capability to efficiently approximate the initial stiff model (2). More precisely, an important step in this process is to calculate the explicit formulae of the ODEs to be solved for the zero-order and the first-order models. This is done by symbolic computations, using Maple [1]. Once the right-hand side terms of the equations are identified, an explicit Runge-Kutta method is implemented, using Maple as a software code generator, in a standard, computationally efficient way to solve the reduced models. The obtained results show that the numerical solution of the first-order two-scale model brings a significant gain in terms of computational time compared to the numerical solution of the full model (2). The reason is that the reduced models use large time steps with respect to the cyclotron period. Finally, we assess the importance of the first-order model over the zero-order one. To give an overview, we illustrate in fig. 1 that the former is more accurate than the latter when compared to a reference solution of equation (2), for two particular instances of the electromagnetic terms. We can thus notice the potential of the first-order two-scale model to follow the drift of the trajectory from the main macroscopic motion in long times, unlike the zero-order model.

The remainder of the paper is as follows. In Section 2 we recall from [11] the general setting for the zero-order and first-order two-scale models approximating a generic singularly perturbed ODE. In Section 3 we detail the formulae for the equation in (2) and we present the strategy for deriving explicit forms by symbolic computation. In Section 4 we give mainly four test cases for which we are able to calculate the explicit forms of the two-scale ODEs. Eventually, in Section 5 we develop numerical simulations to support the approximation results and the efficiency in terms of computational times. More precisely, we demonstrate on several numerical examples that the first-order two-scale solution (i) is much less costly than the reference solution of the initial model and (ii) is much more accurate than the zero-order two-scale solution.

**Notation.** Throughout the paper, {$e_x, e_y, e_z$} stands for the canonical frame of $\mathbb{R}^3$, $I_3$ denotes the $3 \times 3$ identity matrix and $O_3$ the $3 \times 3$ zero matrix. For any vector $v \in \mathbb{R}^3$ we denote its components by $v = (v_x, v_y, v_z)$. For any matrix $A$ we denote its transpose by $A^T$. We consider also the matrices defined for any $\theta \in [0, 2\pi]$ by

$$ R(\theta) = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & \sin \theta \\ 0 & -\sin \theta & \cos \theta \end{pmatrix} \quad \text{and} \quad \mathcal{R}(\theta) = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \sin \theta & 1 - \cos \theta \\ 0 & \cos \theta - 1 & \sin \theta \end{pmatrix}. \quad (3) $$

## 2 A general two-scale framework for stiff ODEs

The aim of this section is to briefly recall the general setting of the two-scale approximation of the general singularly perturbed equation

$$ \frac{d\mathcal{X}_\varepsilon}{dt} = \mathbf{a}(t, \mathcal{X}_\varepsilon) + \frac{1}{\varepsilon}\mathbf{b}(t, \mathcal{X}_\varepsilon), \quad \mathcal{X}_\varepsilon(s) = \mathcal{X}, \quad (4) $$

where $\mathcal{X}_\varepsilon: \mathbb{R} \to \mathbb{R}^6$, $\mathbf{a}$ and $\mathbf{b}$ are given fields satisfying suitable assumptions, and $\mathcal{X} \in \mathbb{R}^6$ is the initial condition. We consider the framework of the two-scale asymptotic expansions to approximate
<figure><img src="image_2.png" /><figcaption>Figure 1: Two trajectories of a particle in the position space, subject to two different inhomogeneous strong magnetic fields. The number $P = \omega_c^{-1}$ denotes the cyclotron period, of order $2\pi\epsilon$. The time step of the zero-order and first-order models is $\Delta t = 0.58 \approx 9\omega_c^{-1}$. The starting position is (1, 1, 1).</figcaption></figure>

the solution $\mathcal{X}_{\epsilon}(t)$ when $\epsilon \to 0$. Specifically, under some regularity assumptions on $\mathbf{a}$ and $\mathbf{b}$ and the $2\pi$-periodicity assumption in the variable $\theta$ on the solution $\mathbf{Z}(t, \theta, \mathbf{z})$ of the equation

$$ \frac{d\mathbf{Z}}{d\theta} = \mathbf{b}(t, \mathbf{Z}), \quad \mathbf{Z}(t, 0, \mathbf{z}) = \mathbf{z}, \qquad (5) $$

for every $\mathbf{z} \in \mathbb{R}^6$ and for every $t \in \mathbb{R}_+$, it has been shown in [11] that $\mathcal{X}_{\epsilon}$ has the following two-scale expansion in time

$$ \mathcal{X}_{\epsilon}(t) = \mathcal{X}^{0}\left(t, \frac{t-s}{\epsilon}\right) + \epsilon \mathcal{X}^{1}\left(t, \frac{t-s}{\epsilon}\right) + \epsilon^{2} \mathcal{X}^{2}\left(t, \frac{t-s}{\epsilon}\right) + \dots \qquad (6) $$

when $\epsilon \to 0$ and where the functions $\mathcal{X}^i(t, \theta)$ are periodic in $\theta$ for every $i \in \mathbb{N}$. Rigorous results on the existence of such expansions are reported in the literature (see, e.g., [24]). The next step is to show that the terms appearing in the expansion (6) are solutions to equations that need to be identified. This was achieved in [11, Theorems 1.1–1.3] for the framework entailed in equation (4). Finally, the last step is to justify the approximations these terms give to the original solution $\mathcal{X}_{\epsilon}$, which is obtained through strong convergence theorems (see [11, Theorems 1.1–1.3]). In the following, we summarise these ideas and recall the expressions of the zero and the first-order approximations, together with the associated convergence results, in the case of equation (4). With
similar notations to [11], we define for any $(Y, Y^0, Y^1) \in (\mathbb{R}^6)^3$ and any $\theta \in [0, 2\pi]$ the following fields

$$ \tilde{\alpha}^0(t, \theta, Y) := \{\nabla_y Z(t, \theta, Y)\}^{-1} \{a(t, Z(t, \theta, Y)) - \partial_t Z(t, \theta, Y)\}, $$

$$ \tilde{\mathbf{a}}^0(t, Y) := \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^0(t, \theta, Y) d\theta, $$

$$ \tilde{\mathbf{A}}^0(t, \theta, Y) := \frac{1}{\theta} \int_0^\theta \tilde{\alpha}^0(t, \sigma, Y) d\sigma - \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^0(t, \sigma, Y) d\sigma, $$

$$ \mathcal{X}^1(t, \theta) := \nabla_y Z(t, \theta, Y^0)(Y^1 + \theta \tilde{\mathbf{A}}^0(t, \theta, Y^0)), $$

$$ \tilde{\alpha}^1(t, \theta, Y^0, Y^1) := \{\nabla_y Z(t, \theta, Y^0)\}^{-1} \left\{ \frac{1}{2} \nabla_y^2 \mathbf{b}(t, Z(t, \theta, Y^0)) \{\mathcal{X}^1(t, \theta); \mathcal{X}^1(t, \theta)\} \right. \\ \qquad \left. + \nabla_y \mathbf{a}(t, Z(t, \theta, Y^0)) \mathcal{X}^1(t, \theta) - \partial_t \nabla_y Z(t, \theta, Y^0) (Y^1 + \theta \tilde{\mathbf{A}}^0(t, \theta, Y^0)) \right. \\ \qquad \left. - \nabla_y^2 Z(t, \theta, Y^0) \{\tilde{\mathbf{a}}^0(t, Y^0); Y^1 + \theta \tilde{\mathbf{A}}^0(t, \theta, Y^0)\} \right\} \\ \qquad - \theta \partial_t \tilde{\mathbf{A}}^0(t, \theta, Y^0) - \theta \nabla_y \tilde{\mathbf{A}}^0(t, \theta, Y^0) \tilde{\mathbf{a}}^0(t, Y^0), $$

$$ \tilde{\mathbf{a}}^1(t, Y^0, Y^1) := \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^1(t, \theta, Y^0, Y^1) d\theta. $$

In the previous expressions, $\nabla_y \mathbf{X}$ stands for the Jacobian matrix of $\mathbf{X}$, while the vector $\nabla_y^2 \mathbf{X}$ is defined by its components

$$ \forall i \in \{1, \dots, 6\}, \quad (\nabla_y^2 \mathbf{X}(y)\{\alpha; \beta\})_i := \sum_{j,k=1}^6 \frac{\partial^2 \mathbf{X}_i}{\partial y_j \partial y_k}(y) \alpha_j \beta_k \quad \forall (y, \alpha, \beta) \in (\mathbb{R}^6)^3, \quad (8) $$

where $\mathbf{X}: \mathbb{R}^6 \to \mathbb{R}^6$ is a general vector field. In this setting, the first two terms in the expansion (6) are given by

$$ \mathcal{X}^0(t, \theta) = \mathbf{Z}(t, \theta, Y^0(t)), \quad (9) $$

$$ \mathcal{X}^1(t, \theta) = \nabla_y Z(t, \theta, Y^0) \{Y^1(t) + \theta \tilde{\mathbf{A}}^0(t, \theta, Y^0(t))\}, \quad (10) $$

where the slow unknowns $Y^0(t)$ and $Y^1(t)$ are solutions to the equations

$$ \frac{dY^0}{dt} = \tilde{\mathbf{a}}^0(t, Y^0), \quad Y^0(s) = \mathcal{X}, \quad (11) $$

and

$$ \frac{dY^1}{dt} = \tilde{\mathbf{a}}^1(t, Y^0, Y^1), \quad Y^1(s) = 0, \quad (12) $$

respectively. More precisely, by [11, Theorems 1.1–1.3], for any $\mathcal{X} \in \mathbb{R}^6$, $s \in \mathbb{R}$, $T \in \mathbb{R}$ and any $\varepsilon > 0$, the solution $\mathcal{X}_\varepsilon(\cdot)$ of (4) exists on $[s, s+T]$, is unique, and is approximated at the zero-order by $\mathcal{X}^0$ given in (9)-(11), in the following sense

$$ \lim_{\varepsilon \to 0} \sup_{t \in [s, s+T]} \left| \mathcal{X}_\varepsilon(t) - \mathcal{X}^0\left(t, \frac{t-s}{\varepsilon}\right) \right| = 0, $$
where $|\cdot|$ is the Euclidean norm on $\mathbb{R}^6$. In addition, if we consider $\mathcal{X}^1$ given by (10)–(12), we have the following first-order approximation

$$ \lim_{\epsilon \to 0} \sup_{t \in [s, s+T]} \left| \frac{1}{\epsilon} \left( \mathcal{X}_{\epsilon}(t) - \mathcal{X}^0 \left( t, \frac{t-s}{\epsilon} \right) \right) - \mathcal{X}^1 \left( t, \frac{t-s}{\epsilon} \right) \right| = 0. $$

We call the zero-order (or the limit) two-scale model the system of equations (9)-(11) and the first-order two-scale model the system (9)–(12). From now on, we will simply refer to the first (zero)-order two-scale model as the first (zero)-order model.

### 3 The case of the Newton–Lorentz equation

In this section we explicit the formulae in (7) in the particular case of the Newton–Lorentz equation in (2). This derivation was initially proposed in [11, Section 3.4]. However, we go back to these calculations in order to detail them in the context of our approach to solve numerically equation (2). The methodology is as follows: first, we derive explicit expressions of the right-hand side terms in (11) and (12), by using the symbolic software Maple [1]. Then, we solve these equations by the 4th order explicit Runge–Kutta method (RK4), which is easy to implement once we are able to evaluate the right-hand sides. Finally, we compute the approximations at order 0 and 1, by using (9) and (10).

The first step can be developed as follows. The system of equations in (2) can be written as in (4) by defining the variable $\mathcal{X}_{\epsilon} = (\mathbf{x}_{\epsilon}, \mathbf{v}_{\epsilon}) \in \mathbb{R}^3 \times \mathbb{R}^3$ and the fields $\mathbf{a}$ and $\mathbf{b}$ by

$$ \mathbf{a}(t, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} \mathbf{v} \\ \mathbf{E}(t, \mathbf{x}) + \mathbf{v} \times \mathbf{B}(\mathbf{x}) \end{pmatrix}, \quad \mathbf{b}(t, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} 0 \\ \mathbf{v} \times \mathbf{m}(t, \mathbf{x}) \end{pmatrix}. \quad (13) $$

We recall that $m$ is assumed of constant norm 1: $m_x(t, \mathbf{x})^2 + m_y(t, \mathbf{x})^2 + m_z(t, \mathbf{x})^2 = 1$ for all $(t, \mathbf{x}) \in \mathbb{R}^+ \times \mathbb{R}^3$. In this case, the asymptotic expansion in (6) writes

$$ \begin{aligned} \mathbf{x}_{\epsilon}(t) &= \mathbf{x}^0 \left( t, \frac{t-s}{\epsilon} \right) + \epsilon \mathbf{x}^1 \left( t, \frac{t-s}{\epsilon} \right) + \dots \\ \mathbf{v}_{\epsilon}(t) &= \mathbf{v}^0 \left( t, \frac{t-s}{\epsilon} \right) + \epsilon \mathbf{v}^1 \left( t, \frac{t-s}{\epsilon} \right) + \dots \end{aligned} \quad (14) $$

We next develop the expressions in (7) in terms of $\mathbf{E}$, $\mathbf{B}$, and $m$. To this end, we introduce notations of type $\mathcal{Y} = (\mathcal{Y}_{\mathbf{x}}, \mathcal{Y}_{\mathbf{v}})$ to separate the position components from the velocity components of a 6-dimensional vector $\mathcal{Y}$. More precisely, we will denote $\mathcal{Y}^0$ by $(\mathbf{y}^0, \mathbf{u}^0)^T$, $\mathcal{X}^0$ by $(\mathbf{x}^0, \mathbf{v}^0)^T$, etc.

#### 3.1 The zero-order model

In the following, we explicit the zero-order terms in the expressions in (7) in the case of $\mathbf{a}$ and $\mathbf{b}$ given above. First, solving the corresponding system (5) leads to

$$ \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} \mathbf{x} \\ A\mathbf{v} \end{pmatrix}, \quad (15) $$
where the matrix $A \equiv A(t, \theta, \mathbf{x})$ is given by

$$A = \begin{pmatrix} (m_y^2 + m_z^2) \cos(\theta) + m_x^2 & m_z \sin(\theta) + m_x m_y (1 - \cos(\theta)) & -m_y \sin(\theta) + m_x m_z (1 - \cos(\theta)) \\ -m_z \sin(\theta) + m_x m_y (1 - \cos(\theta)) & (m_x^2 + m_z^2) \cos(\theta) + m_y^2 & m_x \sin(\theta) + m_y m_z (1 - \cos(\theta)) \\ m_y \sin(\theta) + m_x m_z (1 - \cos(\theta)) & -m_x \sin(\theta) + m_y m_z (1 - \cos(\theta)) & (m_x^2 + m_y^2) \cos(\theta) + m_z^2 \end{pmatrix}. \quad (16)$$

**Remark 3.1.** The matrix $A$ is orthogonal, i.e., $AA^T = I_3$.

Then, we compute

$$\nabla_{\mathbf{x},\mathbf{v}} \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} I_3 & O_3 \\ \nabla_{\mathbf{x}}(A\mathbf{v}) & A \end{pmatrix} \quad \text{and} \quad \mathbf{a}(t, \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v})) = \begin{pmatrix} A\mathbf{v} \\ (E(t, \mathbf{x}) + A\mathbf{v}) \times B(\mathbf{x}) \end{pmatrix}$$

and thus

$$(\nabla_{\mathbf{x},\mathbf{v}} \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}))^{-1} \mathbf{a}(t, \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v})) = \begin{pmatrix} A\mathbf{v} \\ -A^T \nabla_{\mathbf{x}}(A\mathbf{v})A\mathbf{v} + A^T E(t, \mathbf{x}) + \mathbf{v} \times (A^T B(\mathbf{x})) \end{pmatrix}.$$

We also have

$$(\nabla_{\mathbf{x},\mathbf{v}} \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}))^{-1} \partial_t \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} 0 \\ A^T (\partial_t A) \mathbf{v} \end{pmatrix}$$

and therefore, by adding the last two equalities, we obtain $\tilde{\alpha}^0(t, \theta, \mathbf{x}, \mathbf{v})$. As an aside, we note that the term $\theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{x}, \mathbf{v})$ may be simply calculated by automatically integrating $\tilde{\alpha}^0$ with Maple. To conclude this part, the first term in the expansion (14) is given for all $t \ge s$ and for all $\theta \in [0, 2\pi]$ by

$$\begin{pmatrix} \mathbf{x}^0(t, \theta) \\ \mathbf{v}^0(t, \theta) \end{pmatrix} = \begin{pmatrix} \mathbf{y}^0(t) \\ A(t, \theta, \mathbf{y}^0(t)) \mathbf{u}^0(t) \end{pmatrix}, \quad (17)$$

where $(\mathbf{y}^0, \mathbf{u}^0)$ is solution to the initial-value problem

$$\left\{ \begin{aligned} \frac{dy^0}{dt} &= \tilde{\mathbf{a}}^0(t, \mathbf{y}^0, \mathbf{u}^0)_x \\ \frac{du^0}{dt} &= \tilde{\mathbf{a}}^0(t, \mathbf{y}^0, \mathbf{u}^0)_v \\ \mathbf{y}^0(s) &= \mathbf{x}, \quad \mathbf{u}^0(s) = \mathbf{v}, \end{aligned} \right. \quad (18)$$

and

$$\begin{aligned} \tilde{\mathbf{a}}^0(t, \mathbf{x}, \mathbf{v})_x &= \frac{1}{2\pi} \int_0^{2\pi} A\mathbf{v} \, d\theta, \\ \tilde{\mathbf{a}}^0(t, \mathbf{x}, \mathbf{v})_v &= \frac{1}{2\pi} \int_0^{2\pi} \left\{ -A^T \nabla_{\mathbf{x}}(A\mathbf{v})A\mathbf{v} + A^T E(t, \mathbf{x}) + \mathbf{v} \times (A^T B(\mathbf{x})) - A^T (\partial_t A)\mathbf{v} \right\} d\theta. \end{aligned} \quad (19)$$

The right-hand side terms $\tilde{\mathbf{a}}^0(t, \mathbf{x}, \mathbf{v})_x$ and $\tilde{\mathbf{a}}^0(t, \mathbf{x}, \mathbf{v})_v$ are calculated implementing these formulae in Maple. The basic input blocks of this computation are the vectors $E(t, \mathbf{x})$ and $B(\mathbf{x})$ and the matrices $A, \partial_t A, \nabla_{\mathbf{x}}(A\mathbf{v})$, calculated with automatic differentiation commands.
## 3.2 The first-order model

In this section we first decompose, for practical reasons, the first-order model in (12) in a sum of five terms and then, we explicit these terms for the case of **a** and **b** given by (13). We rewrite (12) as

$$ \frac{dY^1}{dt} = \frac{1}{2\pi} \int_{0}^{2\pi} (T_1 + T_2 - T_3 - T_4 - T_5) d\theta, \quad (20) $$

where

$$ T_1 = \frac{1}{2} \{ \nabla_y Z(t, \theta, Y^0) \}^{-1} \nabla_y^2 b(t, Z(t, \theta, Y^0)) \{ X^1(t, \theta); X^1(t, \theta) \}, \quad (21) $$

$$ T_2 = \{ \nabla_y Z(t, \theta, Y^0) \}^{-1} \nabla_y a(t, Z(t, \theta, Y^0)) X^1(t, \theta), \quad (22) $$

$$ T_3 = \{ \nabla_y Z(t, \theta, Y^0) \}^{-1} \partial_t \nabla_y Z(t, \theta, Y^0) (Y^1 + \theta \tilde{A}^0(t, \theta, Y^0)) \quad (23) $$

$$ T_4 = \{ \nabla_y Z(t, \theta, Y^0) \}^{-1} \nabla_y^2 Z(t, \theta, Y^0) \{ \tilde{a}^0(t, Y^0); Y^1 + \theta \tilde{A}^0(t, \theta, Y^0) \} \quad (24) $$

and

$$ T_5 = \theta \partial_t \tilde{A}^0(t, \theta, Y^0) + \theta \nabla_y \tilde{A}^0(t, \theta, Y^0) \tilde{a}^0(t, Y^0). \quad (25) $$

We now develop the expressions above for **a** and **b** given by (13). First, in order to explicit the term $T_1$, we can show by using (8) that

$$ (\nabla_{x,v}^2 b(t, x, v) \{\alpha; \beta\})_x = 0 \quad \text{and} $$

$$ (\nabla_{x,v}^2 b(t, x, v) \{\alpha; \beta\})_v = v \times \begin{pmatrix} \beta_x^T \nabla_x^2 m_x(t, x) \alpha_x \\ \beta_x^T \nabla_x^2 m_y(t, x) \alpha_x \\ \beta_x^T \nabla_x^2 m_z(t, x) \alpha_x \end{pmatrix} + \alpha_v \times (\beta_x \cdot \nabla_x) m(t, x) \\ + \beta_v \times (\alpha_x \cdot \nabla_x) m(t, x), $$

for any $\alpha = (\alpha_x, \alpha_v) \in \mathbb{R}^3 \times \mathbb{R}^3$ and $\beta = (\beta_x, \beta_v) \in \mathbb{R}^3 \times \mathbb{R}^3$, and where $\nabla_x^2 m_*$ denotes the Hessian matrix of $m_*$ (where * stands for $x$, $y$, or $z$). Thus, the term $T_1$ can easily be obtained from

$$ (\nabla_{x,v}^2 b(t, Z(t, \theta, y^0, u^0)) \{(x^1, v^1); (x^1, v^1)\})_x = 0 \quad \text{and} $$

$$ (\nabla_{x,v}^2 b(t, Z(t, \theta, y^0, u^0)) \{(x^1, v^1); (x^1, v^1)\})_v = A u^0 \times \begin{pmatrix} (x^1)^T (\nabla_x^2 m_x) x^1 \\ (x^1)^T (\nabla_x^2 m_y) x^1 \\ (x^1)^T (\nabla_x^2 m_z) x^1 \end{pmatrix} + 2 v^1 \times (x^1 \cdot \nabla_x) m, \quad (26) $$

where $A \equiv A(t, \theta, y^0)$ and $m \equiv m(t, y^0)$. Second, we have

$$ \nabla_{x,v} a(t, x, v) = \begin{pmatrix} O_3 \\ \nabla_x E(t, x) + \nabla_x (v \times B(x)) \end{pmatrix} \begin{pmatrix} I_3 \\ M_B(x) \end{pmatrix}, \quad \text{where} \quad M_B(x) = \begin{pmatrix} 0 & B_z(x) & -B_y(x) \\ -B_z(x) & 0 & B_x(x) \\ B_y(x) & -B_x(x) & 0 \end{pmatrix} $$

and thus, the computation of the term $T_2$ is based on

$$ \nabla_{x,v} a(t, Z(t, \theta, y^0, u^0)) \begin{pmatrix} x^1 \\ v^1 \end{pmatrix} = \begin{pmatrix} v^1 \\ [\nabla_x E(t, y^0) + \nabla_x (A(t, \theta, y^0)u^0 \times B(y^0))] x^1 + M_B(y^0) v^1 \end{pmatrix}. \quad (27) $$
Third, we explicit the term $\mathbf{T}_3$. For any vector $\mathbf{Y} = (\mathbf{Y}_x, \mathbf{Y}_v)^T$ we have

$$ \partial_t \nabla_{\mathbf{x},\mathbf{v}} \mathbf{Z}(t, \theta, \mathbf{y}^0, \mathbf{u}^0) \mathbf{Y} = \begin{pmatrix} 0 \\ \partial_t (\nabla_{\mathbf{x}}(A(t, \theta, \mathbf{y}^0)\mathbf{u}^0)) \mathbf{Y}_x + \partial_t A(t, \theta, \mathbf{y}^0) \mathbf{Y}_v \end{pmatrix} \quad (28) $$

and then, we take $\mathbf{Y}_x = \mathbf{y}^1 + \theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{y}^0, \mathbf{u}^0)_x$ and $\mathbf{Y}_v = \mathbf{u}^1 + \theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{y}^0, \mathbf{u}^0)_v$. Thus, we obtain $\mathbf{T}_3$, by using the building blocks already calculated, as $\theta \tilde{\mathbf{A}}^0$, $\nabla_{\mathbf{x}}(A\mathbf{u}^0)$, etc. Next, the term $\mathbf{T}_4$ can be calculated by means of the general formula (8), as we did in (26) for the term $\mathbf{T}_1$. Indeed, we have

$$ (\nabla_{\mathbf{x},\mathbf{v}}^2 \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}) \{\alpha; \beta\})_{\mathbf{x}} = 0 \quad \text{and} $$

$$ (\nabla_{\mathbf{x},\mathbf{v}}^2 \mathbf{Z}(t, \theta, \mathbf{x}, \mathbf{v}) \{\alpha; \beta\})_{\mathbf{v}} = \begin{pmatrix} \beta_{\mathbf{x}}^T \nabla_{\mathbf{x}}^2 (A \mathbf{v})_x \alpha_{\mathbf{x}} \\ \beta_{\mathbf{x}}^T \nabla_{\mathbf{x}}^2 (A \mathbf{v})_y \alpha_{\mathbf{x}} \\ \beta_{\mathbf{x}}^T \nabla_{\mathbf{x}}^2 (A \mathbf{v})_z \alpha_{\mathbf{x}} \end{pmatrix} + (\beta_{\mathbf{x}} \cdot \nabla_{\mathbf{x}}) A \alpha_{\mathbf{v}} + (\alpha_{\mathbf{x}} \cdot \nabla_{\mathbf{x}}) A \beta_{\mathbf{v}}, $$

for any $\alpha = (\alpha_x, \alpha_v) \in \mathbb{R}^3 \times \mathbb{R}^3$ and $\beta = (\beta_x, \beta_v) \in \mathbb{R}^3 \times \mathbb{R}^3$, and where $(\alpha_x \cdot \nabla_{\mathbf{x}})A$ denotes the matrix whose entries are $\alpha_x \cdot \nabla_{\mathbf{x}} A_{ij}$ for all $i, j \in \{1, 2, 3\}$. Finally, the computation of $\mathbf{T}_5$ is based on Maple differentiation of already calculated terms.

Lastly, by adding up in (20) all the terms determined above, we conclude that $(\mathbf{y}^0, \mathbf{u}^0)$ is the solution of the initial-value problem (18)-(19) and $(\mathbf{y}^1, \mathbf{u}^1)$ is the solution of the initial-value problem

$$ \left\{ \begin{aligned} \frac{d\mathbf{y}^1}{dt} &= \tilde{\mathbf{a}}^1(t, \mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)_x \\ \frac{d\mathbf{u}^1}{dt} &= \tilde{\mathbf{a}}^1(t, \mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)_v \\ \mathbf{y}^1(s) &= 0, \quad \mathbf{u}^1(s) = 0. \end{aligned} \right. \quad (29) $$

Then, the second term (at order $\epsilon$) in the expansion (14) is obtained from $(\mathbf{y}^0(t), \mathbf{u}^0(t))$ and $(\mathbf{y}^1(t), \mathbf{u}^1(t))$ by the formula

$$ \begin{aligned} \begin{pmatrix} \mathbf{x}^1(t, \theta) \\ \mathbf{v}^1(t, \theta) \end{pmatrix} &= \nabla_{\mathbf{x},\mathbf{v}} \mathbf{Z}(t, \theta, \mathbf{y}^0(t), \mathbf{u}^0(t)) \begin{pmatrix} \mathbf{y}^1(t) + \theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{y}^0(t), \mathbf{u}^0(t))_x \\ \mathbf{u}^1(t) + \theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{y}^0(t), \mathbf{u}^0(t))_v \end{pmatrix} \\ &= \begin{pmatrix} \tilde{\mathbf{x}}^1 \\ \nabla_{\mathbf{x}}(A(t, \theta, \mathbf{y}^0(t))\mathbf{u}^0(t))\tilde{\mathbf{x}}^1 + A(t, \theta, \mathbf{y}^0(t))\tilde{\mathbf{v}}^1 \end{pmatrix}, \end{aligned} \quad (30) $$

for all $t \ge s$ and for all $\theta \in [0, 2\pi]$, where we denoted $\tilde{\mathbf{x}}^1 = \mathbf{y}^1(t) + \theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{y}^0(t), \mathbf{u}^0(t))_x$ and $\tilde{\mathbf{v}}^1 = \mathbf{u}^1(t) + \theta \tilde{\mathbf{A}}^0(t, \theta, \mathbf{y}^0(t), \mathbf{u}^0(t))_v$.

# 4 Illustrative examples and their dynamics

We develop in this section several specific examples of charged particle dynamics under a strong magnetic field $\frac{1}{\epsilon}m$ with spatial and temporal inhomogeneities, with increasing difficulty. In each case, we describe the setting of the reference solution, computed with the RK4 scheme with a small
time step. The first example considers an inhomogeneous in the $xOy$ plane and constant in time $m$ term. It was introduced in [11, Section 3.4]. For this setting we analyse two cases for the electric term $\mathbf{E}$. The second example has an additional term in the strong magnetic field, which is slowly varying in time. Finally, the third example is similar to the first, but with spatial inhomogeneities in the three directions. As a preliminary step, we first analyse an example where an explicit solution can be computed, in order to tackle a relevant numerical issue identified in the sequel.

## 4.1 An analytic example

We first consider an example where the solutions to the system (2) and to the first-order model can be written explicitly. As a consequence, we are able to calculate numerical-free relative errors of the first-order approximation and also, the relative error of a numerical approximation of the stiff initial model. More precisely, thanks to the first-order two-scale analytic solution, we can compare its approximation error to that of the RK4 approximation of the stiff model, and thus we can find the time step to be used to achieve the required high accuracy (see Section 5).

Hereafter, the coordinates of the vector $\mathbf{x} \in \mathbb{R}^3$ in the canonical basis $\{\mathbf{e}_x, \mathbf{e}_y, \mathbf{e}_z\}$ are denoted by $(x, y, z)$. We take in (2)

$$m \equiv \mathbf{e}_x, \mathbf{B} \equiv 0, \text{ and } \mathbf{E}(\mathbf{x}) = (-x, 0, -z)^T. \qquad (31)$$

The analytic solution of (2)-(31) is easy to derive. Although the strong magnetic term is constant in space and time, the solution of this system has the interesting feature of possessing three scales in time, a behaviour which is known to be a computing challenge for any classic numerical scheme. Therefore, this example captures the intrinsic difficulties of the oscillatory dynamics of a charged particle in a general inhomogeneous magnetic field. Then, by [11, Theorem 3.3], the first-order approximation is given by

$$\mathbf{G}^1(t) = \begin{pmatrix} \mathbf{x}^0(t, \frac{t-s}{\varepsilon}) \\ \mathbf{v}^0(t, \frac{t-s}{\varepsilon}) \end{pmatrix} + \varepsilon \begin{pmatrix} \mathbf{x}^1(t, \frac{t-s}{\varepsilon}) \\ \mathbf{v}^1(t, \frac{t-s}{\varepsilon}) \end{pmatrix}, \qquad (32)$$

where

$$\begin{pmatrix} \mathbf{x}^0(t, \theta) \\ \mathbf{v}^0(t, \theta) \end{pmatrix} = \begin{pmatrix} \mathbf{y}^0(t) \\ R(\theta)\mathbf{u}^0(t) \end{pmatrix} \quad \text{and} \quad \begin{pmatrix} \mathbf{x}^1(t, \theta) \\ \mathbf{v}^1(t, \theta) \end{pmatrix} = \begin{pmatrix} \mathbf{y}^1(t) + \mathcal{R}(\theta)\mathbf{u}^0(t) \\ R(\theta)\mathbf{u}^1(t) + \mathcal{R}(\theta)\mathbf{E}(\mathbf{y}^0(t)) \end{pmatrix}, \qquad (33)$$

and where $(\mathbf{y}^0(t), \mathbf{u}^0(t))$ and $(\mathbf{y}^1(t), \mathbf{u}^1(t))$ are solutions to initial-value problems that are easy to solve. Specifically, we have

$$\left\{ \begin{aligned} \mathbf{y}^0(t) &= (x \cos(t-s) + \mathbf{v}_x \sin(t-s), y, z)^T, \\ \mathbf{u}^0(t) &= (-x \sin(t-s) + \mathbf{v}_x \cos(t-s), \mathbf{v}_y, \mathbf{v}_z)^T, \end{aligned} \right. \qquad (34)$$

and respectively

$$\left\{ \begin{aligned} \mathbf{y}^1(t) &= (0, -z(t-s), 0)^T, \\ \mathbf{u}^1(t) &= (0, \frac{\mathbf{v}_z}{2}(t-s), -\frac{\mathbf{v}_y}{2}(t-s))^T, \end{aligned} \right. \qquad (35)$$

where $(\mathbf{x}, \mathbf{v})$ is the initial condition in (2).

The relevance of this example for the numerical experiments will be discussed at the beginning of Section 5.
## 4.2 Test case 1

As in [11, Section 3.4], we take in (2) as strong magnetic field the following expression

$$m(\mathbf{x}) = \begin{pmatrix} \frac{-y}{\sqrt{x^2 + y^2}} \\ \frac{x}{\sqrt{x^2 + y^2}} \\ 0 \end{pmatrix}. \qquad (36)$$

We furthermore consider the specific cases of

$$\mathbf{E} \equiv 0 \quad \text{and} \quad \mathbf{B} \equiv \mathbf{e}_z, \qquad (37)$$

and

$$\mathbf{E}(\mathbf{x}) = -\mathbf{x} \quad \text{and} \quad \mathbf{B} \equiv \mathbf{e}_z. \qquad (38)$$

We illustrate these trajectories in the physical space in fig. 2, for $\epsilon = 0.01$, final time $T = 30 \approx 477\omega_c^{-1}$, and initial condition $\mathbf{x} = (1, 1, 1)$, $\mathbf{v} = (1, 0, 0)$. More precisely, the 3d dynamics given by (36)-(37) is plotted in fig. 2a and that given by (36)-(38) is plotted in fig. 2c, together with the previous dynamics. Several cross sections are plotted for both cases in figs. 2d to 2f, 2h and 2i.

**Remark 4.1.** We recall that both zero-order and first-order equations were derived in [11], in this case. Note that we can infer from equations (3.47) in [11, Theorem 3.6] that no motion in the $z$ coordinate is captured by the zero-order model.

**Remark 4.2.** We emphasize that there is a slow drift in the $z$ direction, at the scale $\epsilon$, for the case (37), whereas a drift at the scale 1 occurs for the case defined by (38), as illustrated in fig. 2. Thanks to the previous remark, it is clear that these drifts are missed by the zero-order model. The bottom row of fig. 1 illustrates this statement for the case of (38) and that the first-order approximation accurately captures the drifts.

## 4.3 Test case 2

We now consider the case of a strong magnetic field which is inhomogeneous in space and slowly varying in time. We take in (2)

$$m(t, \mathbf{x}) = \begin{pmatrix} \frac{-y}{\sqrt{x^2 + y^2}} \cos(t/2) \\ \frac{x}{\sqrt{x^2 + y^2}} \cos(t/2) \\ \sin(t/2) \end{pmatrix} \qquad (39)$$

and we fix $\mathbf{E}$ and $\mathbf{B}$ given by (37). The trajectory in the physical space is plotted in figs. 2b and 2e to 2g by taking the same parameters as above.

**Remark 4.3.** We can observe in fig. 2e a slow drift in the $xOy$ plane in the radial direction, of order $\epsilon$. We plot also in fig. 2g the banana orbit, obtained by the projection of the trajectory in the $(x + y, z)$ space. The numerical results underline that, unlike the first-order model, the zero-order one does not capture the slow evolution in time of the banana orbit or, equivalently, of the radius in the $xOy$ plane, but it does capture the averaged macroscopic motion in the $\mathbb{R}^3$ position space.
**Remark 4.4.** To fix the ideas, we derive in Appendix the formulae obtained for the zero-order model, starting from the general equations (18)-(19).

## 4.4 Test case 3

Now, starting from the case in (36) we consider an inhomogeneous strong magnetic field, with a non-zero component in the z direction, given by

$$m(\mathbf{x}) = \begin{pmatrix} \frac{-y}{\sqrt{x^2 + y^2}} \cos \left( \delta \frac{y}{\sqrt{x^2 + y^2}} \right) \\ \frac{x}{\sqrt{x^2 + y^2}} \cos \left( \delta \frac{y}{\sqrt{x^2 + y^2}} \right) \\ \sin \left( \delta \frac{y}{\sqrt{x^2 + y^2}} \right) \end{pmatrix}, \quad (40)$$

where $\delta > 0$ is a small parameter. We take **E** and **B** given by (37). The 3d trajectory in the physical space is plotted in fig. 2b and two cross sections are plotted in figs. 2h and 2i. The features of the first-order and zero-order models, similar to those described in remarks 4.2 and 4.3, have been illustrated in the top row in fig. 1.

**Remark 4.5.** The field in (40) was obtained in every point $(x, y)$ by a rotation in the z direction of the vector $m$ given in (36), of angle $\theta$ around the axis of the position vector, where $\theta$ is a small fraction of the sine of the polar angle of the position vector, namely,

$$\theta = \delta y / \sqrt{x^2 + y^2}, \quad \text{with} \quad \delta = 1/5.$$

# 5 Numerical results

In this section we assess through numerical experiments the accuracy and the computational cost of the first-order models for the four instances of equation (2), introduced in Section 4. With the aim of finding a trade-off between accuracy and low computational cost, we deduce the gain of the first-order models compared to the reference solutions of the original model. Comparison with zero-order models is also performed. The numerical experiments were implemented in a Fortran code and were run on a computer with an Intel Core i7-7600U @ 2.80 GHz processor.

First, in order to find the numerical parameter (the time step) defining the reference solution, we use the example in Section 4.1 to evaluate the accuracy of the RK4 solver of (2)-(31) at a fixed final time, when $\epsilon$ is vanishing. Precisely, we need to investigate how small the time step of the RK4 scheme must be in order to achieve an accuracy which is lower than that of the first-order model. We report in Table 1 the $L^\infty$ norm in time over [0, 2] (a short time scale) and over [0, 10] (a long time scale) of the corresponding errors, for different values of $\epsilon$. The value $N_f$ is the number of the quadrature points for solving the smallest time scale of the cyclotron period. Hereafter we denote this period by $P = 2\pi\epsilon$. For instance, we observe that when $\epsilon = 10^{-4}$, 1560 points per $P$ are required to solve (2), in order to reach at the final time $T = 10$ an error which is smaller than that of the first-order model. A smaller value for $N_f$ would lead to a larger error of the RK4 solver
<figure><img src="image_3.png" /><figcaption>Figure 2: Reference solutions of equation (2) with the fields *m*, **B**, **E** given by (36)-(37) (magenta), by (36)-(38) (blue), by (39)-(37) (orange), and by (40)-(37) (green).</figcaption></figure>
<table><thead><tr><th>ε</th><th>10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>10<sup>-2</sup></th></tr></thead><tbody><tr><td>Number N<sub>f</sub></td><td>1700</td><td>300</td><td>60</td></tr><tr><td>Error(RK4 solver)</td><td>1.07 · 10<sup>-8</sup></td><td>1.10 · 10<sup>-6</sup></td><td>7.01 · 10<sup>-5</sup></td></tr><tr><td>Error(1st-order model)</td><td>1.11 · 10<sup>-8</sup></td><td>1.4 · 10<sup>-6</sup></td><td>9.53 · 10<sup>-5</sup></td></tr></tbody></table>

<table><caption>Table 1: Relative errors at final times $T = 2$ (top table) and $T = 10$ (bottom table) when solving (2)-(31). The initial condition is $(1, 1, 1, 1, 1, 1)$.</caption><thead><tr><th>ε</th><th>10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>10<sup>-2</sup></th></tr></thead><tbody><tr><td>Number N<sub>f</sub></td><td>1560</td><td>290</td><td>50</td></tr><tr><td>Error(RK4 solver)</td><td>7.72 · 10<sup>-8</sup></td><td>6.48 · 10<sup>-6</sup></td><td>7.557 · 10<sup>-4</sup></td></tr><tr><td>Error(1st-order model)</td><td>7.78 · 10<sup>-8</sup></td><td>7.27 · 10<sup>-6</sup></td><td>7.961 · 10<sup>-4</sup></td></tr></tbody></table>

and thus to the inability to estimate the error of the first-order model. Therefore, when applying the RK4 scheme on the time interval $[0, 10]$ to the four previous examples, we choose the following values for the time step: $\delta t = P/60$, $\delta t = P/300$, and $\delta t = P/1600$ when $\varepsilon$ is $\varepsilon = 10^{-2}$, $\varepsilon = 10^{-3}$, and $\varepsilon = 10^{-4}$, respectively. In other words, these choices for $\delta t$ can be justified as follows: if we would keep for all $\varepsilon$ the same number $N_f$ for solving the scale $2\pi\varepsilon$, then the accuracy of the RK4 solver at a fixed final time would decrease as $\varepsilon$ decreases, the reason being that the fixed final time leads to an increasing accumulation of error when $\varepsilon$ decreases.

**Test case 1:** We consider now the first test case in Section 4.2, given by (36)-(37). In fig. 3 (left panel) we plot the $L^\infty$ norm in time over $[0, 10]$ of the absolute errors of the first-order model with respect to the reference solutions defined above. In fig. 3 (right panel) we represent the product of this error with the computational time of the first-order numerical solution until the final simulation time $T = 10$. Smaller is this product, more efficient is the numerical approximation. We infer from fig. 3 that a good trade-off between high accuracy and reduced computational cost is obtained for every $\varepsilon$ by taking $\Delta t = 0.58$.

Next, to detail the behaviour of the relative error, we plot its evolution in time with this fixed time step. First, in the top left panel of fig. 4, we illustrate the basic property of smaller errors with smaller $\varepsilon$ for the first-order solution. Then, in the other panels of fig. 4, we compare the first-order and the zero-order solutions in terms of accuracy. We obtain, as expected, smaller errors of the first-order numerical solution and a smaller gap between these errors when $\varepsilon$ decreases, since the slow drift in the $z$-coordinate is of order $\varepsilon$.

Finally, we investigate the efficiency of the first-order numerical scheme in terms of computational time, see Table 2. We first recall how the fine time step is required to decrease with decreasing $\varepsilon$, to guarantee a sufficiently low error, whereas the coarse time step is fixed. Thus, the run time of the reduced model is constant for different $\varepsilon$, of order $10^{-4}$ seconds for this test simulation.

Then, we consider the test case given by (36)-(38) and we plot in fig. 5 the evolution in time of the errors of the solutions of both reduced models with the fixed time step $\Delta t = 0.58$. This figure is to be compared to fig. 4. First, the first-order numerical solution still satisfies the property of smaller errors with smaller values of $\varepsilon$. Then, we underline the large gap (a factor of 10) between the errors of the zero-order approximations and those of the first-order approximations, showing
<figure><img src="image_4.png" /><figcaption>Figure 3: Comparison between the first-order model computed with several time steps and the stiff equation with several values for $\epsilon$ – 1st test case in Section 4.2. Left panel: Absolute errors of the first-order model with respect to a reference solution of the stiff equation. Right panel: The product of the absolute error by the run time of the first-order model. The simulation time is $T = 10$.</figcaption></figure>

that the latter are much more accurate than the former for this test; the reason was illustrated in the bottom row in fig. 1. We conclude this part by highlighting the robustness in terms of errors of the first-order solutions with respect to the choice of the electric field term $E$.

**Test case 2:** Following the same procedure, we solve in this part the equation (2) with the terms given in Section 4.3. We fix the same numerical parameters as before and we treat the three questions in the same order. First, we choose $\Delta t = 0.58$ to obtain a reduced computational time and enough accuracy of the first-order numerical solution (see fig. 6). Second, we detail in fig. 7 the evolution in time of the relative errors, for three decreasing values of $\epsilon$. We obtain a better accuracy of the first-order model than that of the zero-order model, as expected. Note that the error of the first-order solution is larger than the ones obtained in all the other test cases, most likely due to the complexity of the dynamics. Third, we assess in Table 3 the gain in run time of the first-order model with respect to the reference solution. Similarly as before, the run time of the reduced model is constant for different $\epsilon$, of order $5 \cdot 10^{-3}$ seconds, leading to significant gain factors with respect to the reference solution, for example a factor of order $10^4$ in the case $\epsilon = 10^{-4}$.

**Test case 3:** The same analysis is carried out for the example in Section 4.4, in figs. 8 and 9 and Table 4. We observe the higher accuracy of the first-order numerical solution when compared to that of the zero-order one. In addition, we underline that the magnitude of the errors of the first-order solution is very close to that obtained in the test cases defined in Section 4.2.

Eventually, we remark the increasing computational time of the first-order numerical solution with the ordering of the examples, which is caused by the increasing length and complexity of the right-hand side terms of the ODEs. Nevertheless, the first-order model remains effective when compared to the original model even in the most complex case of Section 4.4 (see Table 4).
<figure><img src="image_5.png" /><figcaption>Figure 4: The evolution in time of the relative errors with respect to the reference solution of the first-order model against those of the zero-order model – 1st test case in Section 4.2. The time step of the reduced models is $\Delta t = 0.58$ and those of the reference solutions are calculated following the first line in Table 2. The top left panel plots the errors of the first-order model for different $\epsilon$.</figcaption></figure>

<table><caption>Table 2: The computation time numbers of the reduced and several fine integrators when solving (2)-(36)-(37) up to the final time $T = 10$. The time step of the 1st-order integrator is $\Delta t = 0.58$.</caption><thead><tr><th>ε</th><th>10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>10<sup>-2</sup></th></tr></thead><tbody><tr><td>Number N<sub>f</sub></td><td>1600</td><td>300</td><td>60</td></tr><tr><td>Time(RK4 solver) in sec.</td><td>23.02</td><td>0.41</td><td>8.5 · 10<sup>-3</sup></td></tr><tr><td>Time(1st-order solver) in sec.</td><td>1.4 · 10<sup>-4</sup></td><td>7.7 · 10<sup>-5</sup></td><td>8.3 · 10<sup>-5</sup></td></tr></tbody></table>

<table><caption>Table 3: The computation time numbers of the reduced and several fine integrators when solving (2)-(39)-(37) up to the final time $T = 10$. The time step of the 1st-order integrator is $\Delta t = 0.58$.</caption><thead><tr><th>ε</th><th>10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>10<sup>-2</sup></th></tr></thead><tbody><tr><td>Number N<sub>f</sub></td><td>1600</td><td>300</td><td>60</td></tr><tr><td>Time(RK4 solver) in sec.</td><td>25.9</td><td>0.50</td><td>1.0 · 10<sup>-2</sup></td></tr><tr><td>Time(1st-order solver) in sec.</td><td>5.1 · 10<sup>-3</sup></td><td>5.2 · 10<sup>-3</sup></td><td>4.6 · 10<sup>-3</sup></td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Figure 5: The evolution in time of the relative errors with respect to the reference solution of the first-order model against those of the zero-order model – 2nd test case in Section 4.2. The time step of the reduced models is $\Delta t = 0.58$ and those of the reference solutions are calculated following the first line in Table 2. The top left panel plots the errors of the first-order model for different $\epsilon$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison between the first-order model computed with several time steps and the stiff equation with several values for $\epsilon$ – test case in Section 4.3. Left panel: Absolute errors of the first-order model with respect to a reference solution of the stiff equation. Right panel: The product of the absolute error by the run time of the first-order model. The simulation time is $T = 10$.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

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

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

<figure><img src="image_11.png" /><figcaption>Figure 7: The evolution in time of the relative errors with respect to the reference solution of the first-order model against those of the zero-order model – test case in Section 4.3. The time step of the reduced models is $\Delta t = 0.58$ and those of the reference solutions are calculated following the first line in Table 3. The top left panel plots the errors of the first-order model for different $\epsilon$.</figcaption></figure>

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

<figure><img src="image_13.png" /><figcaption>Figure 8: Comparison between the first-order model computed with several time steps and the stiff equation with several values for $\epsilon$ – test case in Section 4.4. Left panel: Absolute errors of the first-order model with respect to a reference solution of the stiff equation. Right panel: The product of the absolute error by the run time of the first-order model. The simulation time is $T = 10$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 9: The evolution in time of the relative errors with respect to the reference solution of the first-order model against those of the zero-order model – test case in Section 4.4. The time step of the reduced models is $\Delta t = 0.58$ and those of the reference solutions are calculated following the first line in Table 4. The top left panel plots the errors of the first-order model for different $\epsilon$.</figcaption></figure>

<table><caption>Table 4: The computation time numbers of the reduced and several fine integrators when solving (2)-(40)-(37) up to the final time $T = 10$. The time step of the 1st-order integrator is $\Delta t = 0.58$.</caption><thead><tr><th>ε</th><th>10<sup>-4</sup></th><th>10<sup>-3</sup></th><th>10<sup>-2</sup></th></tr></thead><tbody><tr><td>Number N<sub>f</sub></td><td>1600</td><td>300</td><td>60</td></tr><tr><td>Time(RK4 solver) in sec.</td><td>25.9</td><td>0.47</td><td>9.9 · 10<sup>-3</sup></td></tr><tr><td>Time(1st-order solver) in sec.</td><td>1.42 · 10<sup>-2</sup></td><td>1.30 · 10<sup>-2</sup></td><td>1.39 · 10<sup>-2</sup></td></tr></tbody></table>
# 6 Conclusion

The main goal of the current study was to build an efficient and robust numerical frame for the development of first-order two-scale integrators to solve multi-scale ODEs relevant to plasma physics. We have succeeded to overcome the implementation bottleneck of first-order two-scale models initially proposed in [11]. The numerical experiments underlined the gain in computational time that such a model brings in the approximation of the original stiff solution. Furthermore, our investigation highlighted that the first-order numerical solution, in contrast to the zero-order one, is accurate when following the oscillatory dynamics of a charged particle in a general electromagnetic field. As a further step, we propose to improve the numerical solver, namely the RK4 scheme, to take into account the conservation properties of the dynamics, while preserving the gain in calculation time.

# Appendix

Using Maple, we derive explicit formulae of the right-hand side terms in (18)-(19), in the particular case of the magnetic term $m$ given by (39). First, we have

$$ \tilde{\mathbf{a}}^0(t, \mathbf{x}, \mathbf{v})_{\mathbf{x}} = \frac{1}{2\pi} \int_0^{2\pi} A \mathbf{v} d\theta = \left( \frac{\frac{y(y\mathbf{v}_x - x\mathbf{v}_y)}{\Omega^2} \cos^2(t/2) + \frac{y\mathbf{v}_z}{\Omega} \sin(t)/2}{\frac{x(x\mathbf{v}_y - y\mathbf{v}_x)}{\Omega^2} \cos^2(t/2) + \frac{x\mathbf{v}_z}{\Omega} \sin(t)/2} \right) \frac{x\mathbf{v}_y - y\mathbf{v}_x}{\Omega} \sin(t)/2 + \mathbf{v}_z \sin^2(t/2) \quad (41) $$

where we denoted by $\Omega = \sqrt{x^2 + y^2}$. Second, we decompose

$$ \tilde{\mathbf{a}}^0(t, \mathbf{x}, \mathbf{v})_{\mathbf{v}} = \mathcal{J}^{(1)} + \mathcal{J}^{(2)} + \mathcal{J}^{(3)} + \mathcal{J}^{(4)}, \quad (42) $$

where

$$ \mathcal{J}^{(1)} = -\frac{1}{2\pi} \int_{0}^{2\pi} A^{T} \nabla_{\mathbf{x}}(A\mathbf{v})A\mathbf{v} \, d\theta, $$

$$ \mathcal{J}^{(2)} = \frac{1}{2\pi} \int_{0}^{2\pi} A^{T} \mathbf{E}(t, \mathbf{x}) \, d\theta, $$

$$ \mathcal{J}^{(3)} = \frac{1}{2\pi} \int_{0}^{2\pi} \mathbf{v} \times (A^{T} \mathbf{B}(\mathbf{x})) \, d\theta, $$

$$ \mathcal{J}^{(4)} = -\frac{1}{2\pi} \int_{0}^{2\pi} A^{T}(\partial_{t}A)\mathbf{v} \, d\theta. $$
Then, we compute the components of $\mathcal{J}^{(1)}$ and $\mathcal{J}^{(4)}$

$$ \mathcal{J}^{(1)}_x = -\frac{3}{2\Omega^2}\cos^2(t/2)\left[\mathbf{v}_y(x\mathbf{v}_y - y\mathbf{v}_x)\left(\cos^2(t/2) - \frac{1}{3}\right) + x\mathbf{v}_z^2\sin^2(t/2)\right] \\ -\frac{3}{2\Omega^3}\sin(t)\mathbf{v}_z\left[\left(\mathbf{v}_y(x^2 + \frac{y^2}{2}) - \frac{xy\mathbf{v}_x}{2}\right)\cos^2(t/2) - \frac{\mathbf{v}_y}{6}\Omega^2\right] \quad (43) $$

$$ \mathcal{J}^{(1)}_y = \frac{3}{2\Omega^2}\cos^2(t/2)\left[\mathbf{v}_x(x\mathbf{v}_y - y\mathbf{v}_x)\left(\cos^2(t/2) - \frac{1}{3}\right) - y\mathbf{v}_z^2\sin^2(t/2)\right] \\ + \frac{3}{2\Omega^3}\sin(t)\mathbf{v}_z\left[\left(\mathbf{v}_x(y^2 + \frac{x^2}{2}) - \frac{xy\mathbf{v}_y}{2}\right)\cos^2(t/2) - \frac{\mathbf{v}_x}{6}\Omega^2\right] \quad (44) $$

$$ \mathcal{J}^{(1)}_z = \frac{3}{2\Omega^3}\sin(t/2)\cos^2(t/2)(x\mathbf{v}_x + y\mathbf{v}_y)\left[\mathbf{v}_z\Omega\sin(t/2) + (x\mathbf{v}_y - y\mathbf{v}_x)\cos(t/2)\right] \quad (45) $$

$$ \mathcal{J}^{(4)} = \left( \frac{y\mathbf{v}_z}{2\Omega}, \quad -\frac{x\mathbf{v}_z}{2\Omega}, \quad \frac{x\mathbf{v}_y - y\mathbf{v}_x}{2\Omega} \right)^T \quad (46) $$

Next, we have $\mathcal{J}^{(2)} = m^{(2)}\mathbf{E}(t, \mathbf{x})$ and $\mathcal{J}^{(3)} = m^{(3)}\mathbf{B}(\mathbf{x})$, where the matrices $m^{(2)}$ and $m^{(3)}$ are given by

$$ m^{(2)} = \begin{pmatrix} \frac{y^2}{\Omega^2} \cos^2(t/2) & -\frac{xy}{\Omega^2} \cos^2(t/2) & -\frac{y \sin(t)}{\Omega} \\ -\frac{xy}{\Omega^2} \cos^2(t/2) & \frac{x^2}{\Omega^2} \cos^2(t/2) & \frac{x \sin(t)}{\Omega} \\ -\frac{y \sin(t)}{\Omega} & \frac{x \sin(t)}{\Omega} & \sin^2(t/2) \end{pmatrix}, $$

$$ m^{(3)} = \begin{pmatrix} \frac{xy\mathbf{v}_z}{\Omega^2} \cos^2(t/2) - \frac{y\mathbf{v}_y \sin(t)}{\Omega} & -\frac{x^2\mathbf{v}_z}{\Omega^2} \cos^2(t/2) + \frac{x\mathbf{v}_y \sin(t)}{\Omega} & \mathbf{v}_y \sin^2(t/2) - \frac{x\mathbf{v}_z \sin(t)}{\Omega} \\ \frac{y^2\mathbf{v}_z}{\Omega^2} \cos^2(t/2) + \frac{y\mathbf{v}_x \sin(t)}{\Omega} & -\frac{xy\mathbf{v}_z}{\Omega^2} \cos^2(t/2) - \frac{x\mathbf{v}_x \sin(t)}{\Omega} & -\mathbf{v}_x \sin^2(t/2) - \frac{y\mathbf{v}_z \sin(t)}{\Omega} \\ -\frac{y(x\mathbf{v}_x + y\mathbf{v}_y)}{\Omega^2} \cos^2(t/2) & \frac{x(x\mathbf{v}_x + y\mathbf{v}_y)}{\Omega^2} \cos^2(t/2) & \frac{x\mathbf{v}_x + y\mathbf{v}_y}{\Omega} \sin(t) \end{pmatrix}. $$

Gathering (41) and (42) and replacing in (11) we obtain the equation for the evolution in time of $\mathcal{Y}^0 = (\mathbf{y}^0, \mathbf{u}^0)$. Then, we obtain the zero-order approximation of $(\mathbf{x}_\epsilon(t), \mathbf{v}_\epsilon(t))$ for all $t \ge s$ when $\epsilon \to 0$ by

$$ \begin{pmatrix} \mathbf{x}^0(t, \frac{t-s}{\epsilon}) \\ \mathbf{v}^0(t, \frac{t-s}{\epsilon}) \end{pmatrix} = \begin{pmatrix} \mathbf{y}^0(t) \\ A(t, \frac{t-s}{\epsilon}, \mathbf{y}^0(t)) \mathbf{v}^0(t) \end{pmatrix}, $$

where the matrix $A$ is calculated by replacing $m$ given by (39) in (16).

**Remark 6.1.** It is worth noting that even if these explicit forms seem complicated, they can be obtained with a few lines of Maple code. Then, they are exported in the Fortran code to implement the right-hand side of the ODE to be solved by the RK4 scheme.

## References

[1] Maple (2023). Maplesoft, a division of Waterloo Maple Inc., Waterloo, Ontario.
[2] Y. Barsamian, S. A. Hirstoaga, and E. Violard. Efficient data layouts for a three-dimensional electrostatic particle-in-cell code. *Journal of computational science*, 27:345-356, 2018.

[3] C. K. Birdsall and A. B. Langdon. *Plasma Physics via Computer Simulation*. McGraw-Hill, New York, 1985.

[4] M. Bostan. Gyrokinetic vlasov equation in three dimensional setting. second order approximation. *Multiscale Modeling & Simulation*, 8(5):1923-1957, 2010.

[5] A. J. Brizard and T. S. Hahm. Foundations of nonlinear gyrokinetic theory. *Reviews of modern physics*, 79(2):421-468, 2007.

[6] N. Crouseilles, E. Frenod, S. A. Hirstoaga, and A. Mouton. Two-scale macro-micro decomposition of the vlasov equation with a strong magnetic field. *Mathematical Models and Methods in Applied Sciences*, 23(08):1527-1559, 2013.

[7] V. K. Decyk and T. V. Singh. Particle-in-cell algorithms for emerging computer architectures. *Computer Physics Communications*, 185(3):708-719, 2014.

[8] A. V. Dudkovskia, H. R. Wilson, J. W. Connor, D. Dickinson, and F.I. Parra. Nonlinear second order electromagnetic gyrokinetic theory for a tokamak plasma. *Plasma Physics and Controlled Fusion*, 65(4):045010, 2023.

[9] F. Filbet and L. M. Rodrigues. Asymptotically stable particle-in-cell methods for the vlasov-poisson system with a strong external magnetic field. *SIAM Journal on Numerical Analysis*, 54(2):1120-1146, 2016.

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

[11] E. Frénod. Application of the averaging method to the gyrokinetic plasma. *Asymptotic Analysis*, 46:1-28, 2006.

[12] E. Frénod, M. Gutnic, and S. A. Hirstoaga. First order two-scale particle-in-cell numerical method for the vlasov equation. In *Esaim: Proceedings*, volume 38, pages 348-360. EDP Sciences, 2012.

[13] E. Frénod and M. Lutz. On the geometrical gyro-kinetic theory. *Kinetic Relat. Models*, 7(4):621-659, 2014.

[14] F. Hariri, T.-M. Tran, A. Jocksch, E. Lanti, J. Progsch, P. Messmer, S. Brunner, C. Gheller, and L. Villard. A portable platform for accelerated pic codes and its application to gpus using openacc. *Computer Physics Communications*, 207:69-82, 2016.

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

[16] Y. He, Y. Sun, J. Liu, and H. Qin. Volume-preserving algorithms for charged particle dynamics. *Journal of Computational Physics*, 281:135-147, 2015.

[17] R. W. Hockney and J. W. Eastwood. *Computer Simulation Using Particles*. Institute of Physics, Philadelphia, 1988.
[18] C. Knapp, A. Kendl, A. Koskela, and A. Ostermann. Splitting methods for time integration of trajectories in combined electric and magnetic fields. *Physical Review E*, 92(6):063310, 2015.

[19] W.W Lee. Gyrokinetic particle simulation model. *Journal of Computational Physics*, 72(1):243-269, 1987.

[20] X. Li and B. Wang. Long term analysis of splitting methods for charged-particle dynamics. *Applied Mathematics and Computation*, 441:127682, 2023.

[21] R. G. Littlejohn. A guiding center Hamiltonian: A new approach. *Journal of Mathematical Physics*, 20(12):2445-2458, 1979.

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

[23] L. F. Ricketson and L. Chacón. Asymptotic-preserving gyrokinetic implicit particle-orbit integrator for arbitrary electromagnetic fields. *Journal of Computational Physics*, 513:113136, 2024.

[24] J. A. Sanders, F. Verhulst, and J. Murdock. *Averaging methods in nonlinear dynamical systems*, volume 59. Springer, 2007.

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