# A first-order reduced model for a highly oscillating differential equation with application in Penning traps

Sever Adrian Hirstoaga

► To cite this version:

Sever Adrian Hirstoaga. A first-order reduced model for a highly oscillating differential equation with application in Penning traps. SIAM Journal on Scientific Computing, 2024, pp.S225-S245. 10.1137/23M158351X. http://hal-04361370

HAL Id: hal-04361370

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

Submitted on 22 Dec 2023

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 CC BY 4.0 - Attribution - International License
# A first-order reduced model for a highly oscillating differential equation with application in Penning traps

Sever A. Hirstoaga

December 18, 2023

## Abstract

We derive a reduced first-order model from a two-scale asymptotic expansion in a small parameter, in order to approximate the solution of a stiff differential equation. The problem of interest is a multi-scale Newton-Lorentz equation modeling the dynamics of a charged particle under the influence of a linear electric field and of a perturbed strong magnetic field. First, we show that in short times, the first-order model provides a much better approximation than the zero-order one, since it contains terms evolving at slow time scales. Then, thanks to the source-free property of the equations, we propose a volume-preserving method using a particular splitting technique to solve numerically the first-order model. Finally, it turns out that the first-order model does not systematically provide a satisfactory approximation in long times. To overcome this issue, we implement a recent strategy based on the Parareal algorithm, in which the first-order approximation is used for the coarse solver. This approach allows to perform efficient and accurate long-time simulations for any small parameter. Numerical results for two realistic Penning traps are provided to support these statements.

**Keywords:** multi-scale model, Newton-Lorentz equation, first-order approximation, volume-preserving integrator, Parareal algorithm, speedup

## 1 Introduction

The starting point of the present contribution is an efficient solving of the following six dimensional dynamical system for $0 < \varepsilon \ll 1$

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

where $\mathbf{B}_0$ is a constant magnetic field $\mathbf{B}_0 = \vec{e}_x$, where $\{\vec{e}_x, \vec{e}_y, \vec{e}_z\}$ is the canonical frame of $\mathbb{R}^3$ and $(\mathbf{E}, \mathbf{B}) : \mathbb{R}^3 \to \mathbb{R}^6$ is a given electro-magnetic field. The vector $(\mathbf{x}, \mathbf{v})$ is the initial condition given at the initial time $t = s$. The system in (1) models the dynamics of a charged particle in a given electro-magnetic field with perturbations on the magnetic part. We fix to 1 both the charge and the mass of the particle. In this context, $\mathbf{x}_\varepsilon : \mathbb{R} \to \mathbb{R}^3$ stands for the position unknown and $\mathbf{v}_\varepsilon : \mathbb{R} \to \mathbb{R}^3$ for the velocity unknown. The parameter
$1/\epsilon$ in front of the $B_0$ term means that the constant magnetic field is assumed high with respect to the other terms.

In general, the numerical solving of the stiff system (1) is challenging, since any explicit method needs to solve the smallest time scale, of order $\epsilon$, which can become computationally very demanding (see for instance [7, 3, 17, 10, 9, 14] dealing with similar problems).

In addition, the present work is motivated by a specific application: when the electric field is given by<fnref n="1" />

$$ \mathbf{E}(\mathbf{x}) = c \begin{pmatrix} -\mathbf{x}_x \\ \mathbf{x}_y/2 \\ \mathbf{x}_z/2 \end{pmatrix}, \qquad (2) $$

with an arbitrary constant $c > 0$, the system describes the dynamics of a charged particle in an ideal Penning trap. This is a device for storing charged particles by means of a strong homogeneous magnetic field (the $B_0/\epsilon$ term in (1)) and an inhomogeneous quadrupole electric field [8, 15, 16]. More precisely, the large magnetic field confines a particle radially while the electric field prevents it to go indefinitely and thus providing an axial confinement. However, in a real Penning trap, the magnetic field involves imperfections [15, 16] and therefore, a perturbation term (the role of $\mathbf{B}$ in (1)) has to be added in the equation of motion. In this case, a particle experiences in the axial direction two forces: the one which is due to the electric field and the one due to the $\mathbf{B}$ term. As in [7], in this work, we consider two cases of an inhomogeneous perturbation of the magnetic field: first,

$$ \mathbf{B}(\mathbf{x}) = k \begin{pmatrix} \mathbf{x}_x^2 - \frac{\mathbf{x}_y^2 + \mathbf{x}_z^2}{2} \\ -\mathbf{x}_x \mathbf{x}_y \\ -\mathbf{x}_x \mathbf{x}_z \end{pmatrix}, \qquad (3) $$

corresponding to a magnetic bottle, and second,

$$ \mathbf{B}(\mathbf{x}) = k \begin{pmatrix} \mathbf{x}_z - \mathbf{x}_y \\ \mathbf{x}_z - \mathbf{x}_x \\ \mathbf{x}_x + \mathbf{x}_y \end{pmatrix}, \qquad (4) $$

corresponding to an asymmetric magnetic field, where $k > 0$ is a fixed constant.

In this challenging context, we first derive reduced models approximating the equation (1) when $\epsilon$ vanishes, on the basis of a two-scale asymptotic expansion [2]. Both zero-order and first-order in $\epsilon$ approximations are obtained via this procedure. The advantage of reduced models is, from the numerical point of view, their low computational cost since the time step can be chosen independent of $\epsilon$. Next, we show that the first-order model is much more accurate than the zero-order model, when approximating the initial stiff equation (1). Specifically, recall that a particle in an electromagnetic field like (2)-(3) or (2)-(4) evolves at three time scales [8]: a rapid one corresponding to a rotation around the magnetic field line $B_0$, a slower one associated to a bounce motion parallel to $B_0$, and a much slower scale

<footnote n="1">¹Hereafter, we denote the coordinates of a vector $\mathbf{x} \in \mathbb{R}^3$ in the canonical basis by $(\mathbf{x}_x, \mathbf{x}_y, \mathbf{x}_z)$.</footnote>
corresponding to the electric drift across the field line $B_0$. In this context, we obtain that (i) both reduced models average the fastest rotation motion, (ii) the first-order model is more accurate than the zero-order one in the approximation of the bounce motion, and (iii) the electric drift is completely missed by the zero-order model, unlike the first-order one.

Secondly, we propose a volume-preserving method for the numerical solving of the reduced models. This is a natural choice for achieving at the discrete level the conservation of the invariants of the system [4, 1, 7, 6]. In addition, we develop a specific splitting technique, by which the total evolution vector field of the differential equation is decomposed in several components, each of which is simpler to solve than the initial system. This approach allows us to perform simulations of realistic Penning traps which are accurate in short times. However, as theoretically expected, the reduced model's performance degrades when the final time increases. In order to deal with this issue, we propose to implement a strategy based on the Parareal algorithm [11], that we recently developed in a similar context, in [3]. The strategy is based on using a reduced model to define the coarse solver, since such a model has a low computational cost when compared to the full problem. In addition, we highlight that the Parareal algorithm allows to successfully perform accurate and efficient simulations in long times with the help of the first-order model. If the zero-order approximation is used in the Parareal approach, the efficiency is lost in long runs, since this model is inaccurate in large times. The idea of using a reduced model as coarse solver in the Parareal approach goes back at least to [12]. We also mention the works in [9, 5, 14] that have implemented this method for solving similar singularly perturbed ODEs by using only limit (zero-order) models for the coarse solving. These models were obtained either by classic averaging or by projection onto a slow manifold.

The remainder of the paper is organized as follows. We derive the first-order model in Section 2. The discretization method is presented in Section 3 and the numerical results with the application in Penning traps are developed in Section 4. We first assess in Section 4.1 the accuracy of the reduced models and we measure the savings that these models bring in calculation time. Then, in Section 4.2 we investigate the efficiency of the Parareal algorithm in terms of ideal speedup.

**Notation.** Throughout the paper, $I_3$ denotes the $3 \times 3$ identity matrix, $O_3$ the $3 \times 3$ zero matrix, and $O_n$ is the zero vector of $\mathbb{R}^n$. For any vector $\mathbf{v} = (\mathbf{v}_x, \mathbf{v}_y, \mathbf{v}_z) \in \mathbb{R}^3$ we denote its perpendicular part by $\mathbf{v}_{\perp} = (\mathbf{v}_y, \mathbf{v}_z)$. 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 (5)$$
# 2 Derivation of the first-order reduced model

## 2.1 The general framework of a singularly perturbed dynamical system

In this section we introduce the general singularly perturbed dynamical system, which is a particular case of the system considered in [2],

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

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 given at the initial time $s$. The framework of the two-scale asymptotic expansions can be utilized in this context to approximate the solution $\mathcal{X}_\varepsilon(t)$ when $\varepsilon \to 0$. More precisely, under regularity assumptions on $\mathbf{a}$ and $\mathbf{b}$ and the assumption that the solution $\mathbf{Z}(\theta, \mathbf{z})$ to equation

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

is $2\pi$-periodic for every $\mathbf{z} \in \mathbb{R}^6$, it is proved in [2] that $\mathcal{X}_\varepsilon$ admits the following two-scale expansion in time

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

when $\varepsilon \to 0$ and where the functions $\mathcal{X}^i(t, \theta)$ are periodic in $\theta$ for every $i \in \mathbb{N}$. Such an expansion may be rigourously developed up to any order if enough regularity is imposed on the fields. More precisely, for a given integer $m \ge 1$, the existence of continuous and bounded derivatives of the fields ($\mathbf{a}, \mathbf{b}$) up to the orders ($m, m+1$) allows justifying the approximation given by the first $m$ terms in the expansion in (8). In this setting, ordinary differential equations characterizing the terms of the expansion are derived in [2, Theorems 1.1–1.3]. In addition, strong convergence theorems are proved, justifying the approximation results asserting that, e.g., at the first order,

$$ \mathcal{X}_{\varepsilon}(t) \approx \mathcal{X}^{0}\left(t, \frac{t-s}{\varepsilon}\right) + \varepsilon \mathcal{X}^{1}\left(t, \frac{t-s}{\varepsilon}\right), \quad \text{when } \varepsilon \to 0. $$

Next, we detail these ideas by stating the formulas of the zero and the first order approximations and the convergence results, in the case of equation (6). Using the same notation as in [2], we define for any $(\mathcal{X}, \mathcal{X}^0, \mathcal{X}^1) \in (\mathbb{R}^6)^3$ and any $\theta \in [0, 2\pi]$ the following fields

$$ \begin{aligned} \tilde{\alpha}^0(\theta, \mathcal{X}) &:= \{\nabla_x \mathbf{Z}(\theta, \mathcal{X})\}^{-1} \mathbf{a}(\mathbf{Z}(\theta, \mathcal{X})), \\ \tilde{\mathbf{a}}^0(\mathcal{X}) &:= \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^0(\theta, \mathcal{X}) d\theta, \\ \tilde{\mathbf{A}}^0(\theta, \mathcal{X}) &:= \frac{1}{\theta} \int_0^{\theta} \tilde{\alpha}^0(\sigma, \mathcal{X}) d\sigma - \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^0(\sigma, \mathcal{X}) d\sigma, \\ \tilde{\alpha}^1(\theta, \mathcal{X}^0, \mathcal{X}^1) &:= \{\nabla_x \mathbf{Z}(\theta, \mathcal{X}^0)\}^{-1} \{\nabla \mathbf{a}(\mathbf{Z}(\theta, \mathcal{X}^0)) \nabla_x \mathbf{Z}(\theta, \mathcal{X}^0) \{\mathcal{X}^1 + \theta \tilde{\mathbf{A}}^0(\theta, \mathcal{X}^0)\} \\ &\quad - \theta \nabla_x \tilde{\mathbf{A}}^0(\theta, \mathcal{X}^0) \{\tilde{\mathbf{a}}^0(\mathcal{X}^0)\}, \\ \tilde{\mathbf{a}}^1(\mathcal{X}^0, \mathcal{X}^1) &:= \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^1(\theta, \mathcal{X}^0, \mathcal{X}^1) d\theta. \end{aligned} \qquad (9) $$
In this setting, the first two terms in the expansion (8) are given by

$$X^0(t, \theta) = Z(\theta, Y^0(t)), \qquad (10)$$

$$X^1(t, \theta) = \nabla_x Z(\theta, Y^0)\{Y^1(t) + \theta \tilde{A}^0(\theta, Y^0(t))\}, \qquad (11)$$

where the slow (i.e. independent of $\theta$) quantities $Y^0(t)$ and $Y^1(t)$ are solutions of the equations

$$\frac{dY^0}{dt} = \tilde{a}^0(Y^0), \quad Y^0(s) = X, \qquad (12)$$

and

$$\frac{dY^1}{dt} = \tilde{a}^1(Y^0, Y^1), \quad Y^1(s) = 0, \qquad (13)$$

respectively. More precisely, by [2, Theorems 1.1-1.3], for any $X \in \mathbb{R}^6$, $s \in \mathbb{R}$, $T \in \mathbb{R}$ and any $\varepsilon > 0$, the solution $X_{\varepsilon}(\cdot)$ of (6) exists on $[s, s+T]$, is unique, and is approximate by $X^0$ given in (10)-(12), in the sense of

$$\lim_{\varepsilon \to 0} \sup_{t \in [s, s+T]} \left| X_{\varepsilon}(t) - 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 $X^1$ given by (11)-(13), we have the following first-order approximation

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

## 2.2 The case of Newton-Lorentz equation

The system of equations in (1) is of the form of (6) by taking the variable $X_{\varepsilon} = (x_{\varepsilon}, v_{\varepsilon}) \in (\mathbb{R}^3)^2$ and the fields a and b by

$$a(x, v) = \begin{pmatrix} v \\ E(x) + v \times B(x) \end{pmatrix}, \quad b(v) = \begin{pmatrix} 0 \\ v \times B_0 \end{pmatrix}. \qquad (14)$$

In this case, the asymptotic expansion in (8) writes

$$\begin{aligned} x_{\varepsilon}(t) &= x^0\left(t, \frac{t-s}{\varepsilon}\right) + \varepsilon x^1\left(t, \frac{t-s}{\varepsilon}\right) + \dots \\ v_{\varepsilon}(t) &= v^0\left(t, \frac{t-s}{\varepsilon}\right) + \varepsilon v^1\left(t, \frac{t-s}{\varepsilon}\right) + \dots \end{aligned} \qquad (15)$$

Note that in all the applications discussed in [2, Section 3] the term $v \times B(x)$ in the definition of a in (14) is equal to zero. Therefore, the reduced models previously obtained do not account for this term.

In this work we consider a non-zero $B(x)$ term motivated by the applications to the real Penning traps, see (3)-(4). Thus, in a first stage, the aim is to derive the reduced models of order one from the general formulas in [2, Section 1].

The main theoretical result of this paper is the following.
**Theorem 2.1.** If E and B are in $C_b^2(\mathbb{R}^3)$, the first two terms of the expansion (15) of the solution $(\mathbf{x}_\epsilon(t), \mathbf{v}_\epsilon(t))$ to (1) are 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) \\ R(\theta) \mathbf{u}^0(t) \end{pmatrix} \qquad (16) $$

and

$$ \begin{pmatrix} \mathbf{x}^1(t, \theta) \\ \mathbf{v}^1(t, \theta) \end{pmatrix} = \begin{pmatrix} \mathbf{y}^1(t) \\ R(\theta) \mathbf{u}^1(t) \end{pmatrix} + \begin{pmatrix} \mathcal{R}(\theta) \mathbf{u}^0(t) \\ \mathcal{R}(\theta) \mathbf{E}(\mathbf{y}^0(t)) + R(\theta) \mathbf{u}^0(t) \times \mathcal{R}(\theta) \mathbf{B}(\mathbf{y}^0(t)) \end{pmatrix}, \quad (17) $$

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

$$ \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{y}^0}{dt} = \begin{pmatrix} \mathbf{u}_x^0 \\ 0 \\ 0 \end{pmatrix}, \quad \frac{d\mathbf{u}^0}{dt} = \begin{pmatrix} \mathbf{E}_x(\mathbf{y}^0) \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ \mathbf{B}_x(\mathbf{y}^0) \mathbf{u}_z^0 \\ -\mathbf{B}_x(\mathbf{y}^0) \mathbf{u}_y^0 \end{pmatrix} \\[1em] \mathbf{y}^0(s) = \mathbf{x}, \quad \mathbf{u}^0(s) = \mathbf{v}, \end{array} \right. \qquad (18) $$

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

$$ \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{y}^1}{dt} = \begin{pmatrix} \mathbf{u}_x^1 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ \mathbf{E}_z(\mathbf{y}^0) \\ -\mathbf{E}_y(\mathbf{y}^0) \end{pmatrix} + \begin{pmatrix} \mathbf{u}_y^0 \mathbf{B}_y(\mathbf{y}^0) + \mathbf{u}_z^0 \mathbf{B}_z(\mathbf{y}^0) \\ \mathbf{u}_y^0 \mathbf{B}_x(\mathbf{y}^0) + \mathbf{u}_x^0 \mathbf{B}_y(\mathbf{y}^0) \\ \mathbf{u}_z^0 \mathbf{B}_x(\mathbf{y}^0) + \mathbf{u}_x^0 \mathbf{B}_z(\mathbf{y}^0) \end{pmatrix} \\[1em] \displaystyle \frac{d\mathbf{u}^1}{dt} = \begin{pmatrix} \nabla \mathbf{E}_x(\mathbf{y}^0) \cdot \mathbf{y}^1 \\ 0 \\ 0 \end{pmatrix} + \begin{pmatrix} 0 \\ \mathbf{B}_x(\mathbf{y}^0) \mathbf{u}_z^1 \\ -\mathbf{B}_x(\mathbf{y}^0) \mathbf{u}_y^1 \end{pmatrix} + \begin{pmatrix} 0 \\ \mathbf{u}_z^0 \nabla \mathbf{B}_x(\mathbf{y}^0) \cdot \mathbf{y}^1 \\ -\mathbf{u}_y^0 \nabla \mathbf{B}_x(\mathbf{y}^0) \cdot \mathbf{y}^1 \end{pmatrix} + \mathcal{T}(\mathbf{y}^0, \mathbf{u}^0), \\[1em] \mathbf{y}^1(s) = 0, \quad \mathbf{u}^1(s) = 0, \end{array} \right. \qquad (19) $$

where $\mathcal{T}$ is given by

$$ \mathcal{T}(\mathbf{y}^0, \mathbf{u}^0) = \mathcal{T}_\delta(\mathbf{y}^0, \mathbf{u}^0) + \mathcal{T}_m(\mathbf{y}^0, \mathbf{u}^0) + \begin{pmatrix} \mathbf{B}_y(\mathbf{y}^0)\mathbf{E}_y(\mathbf{y}^0) + \mathbf{B}_z(\mathbf{y}^0)\mathbf{E}_z(\mathbf{y}^0) \\ \mathbf{B}_x(\mathbf{y}^0)\mathbf{E}_y(\mathbf{y}^0) + \mathbf{B}_y(\mathbf{y}^0)\mathbf{E}_x(\mathbf{y}^0) \\ \mathbf{B}_x(\mathbf{y}^0)\mathbf{E}_z(\mathbf{y}^0) + \mathbf{B}_z(\mathbf{y}^0)\mathbf{E}_x(\mathbf{y}^0) \end{pmatrix}, \quad (20) $$
and where the electric and magnetic terms, $\mathcal{T}_e$ and $\mathcal{T}_m$, are respectively given by

$$ \mathcal{T}_e(\mathbf{y}^0, \mathbf{u}^0) = \begin{pmatrix} \mathbf{u}_z^0 \partial_y \mathbf{E}_x(\mathbf{y}^0) - \mathbf{u}_y^0 \partial_z \mathbf{E}_x(\mathbf{y}^0) \\ \mathbf{u}_x^0 \partial_x \mathbf{E}_z(\mathbf{y}^0) \\ -\mathbf{u}_x^0 \partial_x \mathbf{E}_y(\mathbf{y}^0) \end{pmatrix} + \frac{1}{2} (\partial_y \mathbf{E}_y(\mathbf{y}^0) + \partial_z \mathbf{E}_z(\mathbf{y}^0)) \begin{pmatrix} 0 \\ -\mathbf{u}_z^0 \\ \mathbf{u}_y^0 \end{pmatrix} \\ + \frac{1}{2} (\partial_z \mathbf{E}_y(\mathbf{y}^0) - \partial_y \mathbf{E}_z(\mathbf{y}^0)) \begin{pmatrix} 0 \\ \mathbf{u}_y^0 \\ \mathbf{u}_z^0 \end{pmatrix}, $$

$$ \mathcal{T}_{m_x}(\mathbf{y}^0, \mathbf{u}^0) = \mathbf{B}_x(\mathbf{y}^0)\mathbf{B}_z(\mathbf{y}^0)\mathbf{u}_y^0 - \mathbf{B}_x(\mathbf{y}^0)\mathbf{B}_y(\mathbf{y}^0)\mathbf{u}_z^0 \\ -\mathbf{u}_x^0\mathbf{u}_y^0\partial_x\mathbf{B}_y(\mathbf{y}^0) - \mathbf{u}_x^0\mathbf{u}_z^0\partial_x\mathbf{B}_z(\mathbf{y}^0) + \frac{1}{2}(\partial_y\mathbf{B}_y(\mathbf{y}^0) + \partial_z\mathbf{B}_z(\mathbf{y}^0))(\mathbf{u}_y^{02} + \mathbf{u}_z^{02}), $$

$$ \mathcal{T}_{m_y}(\mathbf{y}^0, \mathbf{u}^0) = -\mathbf{B}_x(\mathbf{y}^0)\mathbf{B}_z(\mathbf{y}^0)\mathbf{u}_x^0 + \frac{1}{2}\mathbf{u}_z^0(\mathbf{B}_y(\mathbf{y}^0)^2 + \mathbf{B}_z(\mathbf{y}^0)^2) + (\mathbf{u}_z^0)^2\partial_y\mathbf{B}_x(\mathbf{y}^0) \\ + (\mathbf{u}_x^0)^2\partial_x\mathbf{B}_y(\mathbf{y}^0) - \mathbf{u}_y^0\mathbf{u}_z^0\partial_z\mathbf{B}_x(\mathbf{y}^0) - \frac{1}{2}(\partial_y\mathbf{B}_y(\mathbf{y}^0) + \partial_z\mathbf{B}_z(\mathbf{y}^0))\mathbf{u}_x^0\mathbf{u}_y^0 \\ -\frac{1}{2}(\partial_z\mathbf{B}_y(\mathbf{y}^0) - \partial_y\mathbf{B}_z(\mathbf{y}^0))\mathbf{u}_x^0\mathbf{u}_z^0, $$

$$ \mathcal{T}_{m_z}(\mathbf{y}^0, \mathbf{u}^0) = \mathbf{B}_x(\mathbf{y}^0)\mathbf{B}_y(\mathbf{y}^0)\mathbf{u}_x^0 - \frac{1}{2}\mathbf{u}_y^0(\mathbf{B}_y(\mathbf{y}^0)^2 + \mathbf{B}_z(\mathbf{y}^0)^2) + (\mathbf{u}_y^0)^2\partial_z\mathbf{B}_x(\mathbf{y}^0) \\ + (\mathbf{u}_x^0)^2\partial_x\mathbf{B}_z(\mathbf{y}^0) - \mathbf{u}_y^0\mathbf{u}_z^0\partial_y\mathbf{B}_x(\mathbf{y}^0) - \frac{1}{2}(\partial_y\mathbf{B}_y(\mathbf{y}^0) + \partial_z\mathbf{B}_z(\mathbf{y}^0))\mathbf{u}_x^0\mathbf{u}_z^0 \\ + \frac{1}{2}(\partial_z\mathbf{B}_y(\mathbf{y}^0) - \partial_y\mathbf{B}_z(\mathbf{y}^0))\mathbf{u}_x^0\mathbf{u}_y^0. $$

*Proof:* First, we specify the expressions of the fields in (9) in the particular case of **a** and **b** given by (14). Thus, the solution to the rapid rotation given by (7) is the function

$$ \mathbf{Z}(\theta, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} \mathbf{x} \\ R(\theta)\mathbf{v} \end{pmatrix}, \qquad (21) $$

which is $2\pi$-periodic in $\theta$. In addition, we have for all $(\mathbf{x}, \mathbf{v})$,

$$ \nabla_{(\mathbf{x},\mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}, \mathbf{v}) = \begin{pmatrix} I_3 & O_3 \\ O_3 & R(\theta) \end{pmatrix} \quad \text{and} \quad \{\nabla_{(\mathbf{x},\mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}, \mathbf{v})\}^{-1} = \nabla_{(\mathbf{x},\mathbf{v})} \mathbf{Z}(-\theta, \mathbf{x}, \mathbf{v}). $$

Then, we can easily compute

$$ \tilde{\alpha}^0(\theta, \mathbf{x}, \mathbf{v}) = \{\nabla_{(\mathbf{x},\mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}, \mathbf{v})\}^{-1} \mathbf{a}(\mathbf{Z}(\theta, \mathbf{x}, \mathbf{v})) = \begin{pmatrix} R(\theta)\mathbf{v} \\ R(-\theta)\mathbf{E}(\mathbf{x}) + \mathbf{v} \times R(-\theta)\mathbf{B}(\mathbf{x}) \end{pmatrix}, \\ \tilde{\mathbf{a}}^0(\mathbf{x}, \mathbf{v}) = \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}^0(\theta, \mathbf{x}, \mathbf{v}) d\theta = \begin{pmatrix} \mathbf{v}_{\parallel} \\ \mathbf{E}(\mathbf{x})_{\parallel} + \mathbf{v} \times \mathbf{B}(\mathbf{x})_{\parallel} \end{pmatrix}, \qquad (22) $$

where $\mathbf{v}_{\parallel}$ stands for $\mathbf{v}_{\parallel} = (\mathbf{v} \cdot \vec{e}_x)\vec{e}_x = (\mathbf{v}_x, 0, 0)^T$, and thus

$$ \theta \tilde{\mathbf{A}}^0(\theta, \mathbf{x}, \mathbf{v}) = \left( \int_0^\theta d\sigma - \frac{\theta}{2\pi} \int_0^{2\pi} d\sigma \right) \tilde{\alpha}^0(\sigma, \mathbf{x}, \mathbf{v}) \\ = \left( -\mathcal{R}(\theta)\mathbf{E}(\mathbf{x}) + \mathbf{v} \times (-\mathcal{R}(\theta)\mathbf{B}(\mathbf{x})) \right). \qquad (23) $$
Therefore, at this stage, we can already obtain the equations (16) and (18) from (10) and (12), respectively. Then, starting from (11) we can derive the formula (17) of the first order term $(\mathbf{x}^1, \mathbf{v}^1)$ by using the fact that $R(\theta)(-\mathcal{R}(-\theta)) = \mathcal{R}(\theta)$. We are left with deriving the equations in (19) and (20) from the general equations (13) and (9).

The remainder of the proof consists of two steps: first, we make the field $\tilde{\alpha}^1$ in (9) explicit in our case and second, we average in $\theta$ the result in order to obtain the field $\tilde{\mathbf{a}}^1$. Thus, we start by developing $\tilde{\alpha}^1$ in the case of the fields given in (14). We compute the jacobian of $\mathbf{a}$

$$ \nabla \mathbf{a}(\mathbf{x}, \mathbf{v}) = \begin{pmatrix} O_3 & I_3 \\ \nabla \mathbf{E}(\mathbf{x}) + \nabla_{\mathbf{x}}(\mathbf{v} \times \mathbf{B}(\mathbf{x})) & M_{\mathbf{B}}(\mathbf{x}) \end{pmatrix}, $$

where $\nabla \mathbf{E}(\mathbf{x})$ is the jacobian matrix of the field $\mathbf{E}$ and where the matrix $M_{\mathbf{B}}(\mathbf{x})$ is defined by

$$ M_{\mathbf{B}}(\mathbf{x}) = \begin{pmatrix} 0 & \mathbf{B}_z(\mathbf{x}) & -\mathbf{B}_y(\mathbf{x}) \\ -\mathbf{B}_z(\mathbf{x}) & 0 & \mathbf{B}_x(\mathbf{x}) \\ \mathbf{B}_y(\mathbf{x}) & -\mathbf{B}_x(\mathbf{x}) & 0 \end{pmatrix}, $$

whose utility is in the rewriting of the term $\mathbf{v} \times \mathbf{B}$ as the matrix-vector product $M_{\mathbf{B}}\mathbf{v}$. Therefore, from the computations above, we have

$$ \begin{aligned} & \nabla \mathbf{a}(\mathbf{Z}(\theta, \mathbf{x}^0, \mathbf{v}^0)) \nabla_{(\mathbf{x}, \mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}^0, \mathbf{v}^0) \{ (\mathbf{x}^1, \mathbf{v}^1)^T + \theta \tilde{\mathbf{A}}^0(\theta, \mathbf{x}^0, \mathbf{v}^0) \} \\ & = \begin{pmatrix} O_3 & I_3 \\ \nabla \mathbf{E}(\mathbf{x}^0) + \nabla_{\mathbf{x}}(R(\theta)\mathbf{v}^0 \times \mathbf{B}(\mathbf{x}^0)) & M_{\mathbf{B}}(\mathbf{x}^0) \end{pmatrix} \begin{pmatrix} \mathbf{x}^1 + \mathcal{R}(\theta)\mathbf{v}^0 \\ R(\theta)(\mathbf{v}^1 - \mathcal{R}(-\theta)\mathbf{E}(\mathbf{x}^0) + \mathbf{v}^0 \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0))) \end{pmatrix} \end{aligned} $$

and thus, the first term of $\tilde{\alpha}^1$, denoted by $T^1 = (T_x^1, T_v^1)$, is obtained by applying $\{\nabla_{(\mathbf{x}, \mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}^0, \mathbf{v}^0)\}^{-1}$ to the previous expression:

$$ T^1 := \{\nabla_{(\mathbf{x}, \mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}^0, \mathbf{v}^0)\}^{-1} \{\nabla \mathbf{a}(\mathbf{Z}(\theta, \mathbf{x}^0, \mathbf{v}^0)) \nabla_{(\mathbf{x}, \mathbf{v})} \mathbf{Z}(\theta, \mathbf{x}^0, \mathbf{v}^0) \{ (\mathbf{x}^1, \mathbf{v}^1)^T + \theta \tilde{\mathbf{A}}^0(\theta, \mathbf{x}^0, \mathbf{v}^0) \} \}, $$

and we compute

$$ T_x^1 = R(\theta) \left( \mathbf{v}^1 - \mathcal{R}(-\theta) \mathbf{E}(\mathbf{x}^0) + \mathbf{v}^0 \times (-\mathcal{R}(-\theta) \mathbf{B}(\mathbf{x}^0)) \right) $$

and

$$ T_v^1 = \begin{aligned} & (R(-\theta)\nabla\mathbf{E}(\mathbf{x}^0) + R(-\theta)\nabla_{\mathbf{x}}(R(\theta)\mathbf{v}^0 \times \mathbf{B}(\mathbf{x}^0))) [\mathbf{x}^1 + \mathcal{R}(\theta)\mathbf{v}^0] \\ & + R(-\theta)M_{\mathbf{B}}(\mathbf{x}^0)R(\theta) (\mathbf{v}^1 - \mathcal{R}(-\theta)\mathbf{E}(\mathbf{x}^0) + \mathbf{v}^0 \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0))) \end{aligned} $$

By using

$$ R(\theta)(-\mathcal{R}(-\theta)) = \mathcal{R}(\theta) \quad \text{and} \quad R(\theta)(\mathbf{u} \times \mathbf{v}) = R(\theta)\mathbf{u} \times R(\theta)\mathbf{v}, \quad \forall \mathbf{u}, \mathbf{v} \in \mathbb{R}^3 \quad (24) $$

and the fact that $M_{\mathbf{B}}(\mathbf{x})\mathbf{v} = \mathbf{v} \times \mathbf{B}(\mathbf{x})$ for any $\mathbf{v} \in \mathbb{R}^3$, the term $T^1$ takes the following form

$$ T_x^1 = R(\theta)\mathbf{v}^1 + \mathcal{R}(\theta)\mathbf{E}(\mathbf{x}^0) + R(\theta)\mathbf{v}^0 \times \mathcal{R}(\theta)\mathbf{B}(\mathbf{x}^0), \quad (25) $$

$$ T_v^1 = \begin{aligned} & (R(-\theta)\nabla\mathbf{E}(\mathbf{x}^0) + R(-\theta)\nabla_{\mathbf{x}}(R(\theta)\mathbf{v}^0 \times \mathbf{B}(\mathbf{x}^0))) [\mathbf{x}^1 + \mathcal{R}(\theta)\mathbf{v}^0] \\ & + (\mathbf{v}^1 - \mathcal{R}(-\theta)\mathbf{E}(\mathbf{x}^0) + \mathbf{v}^0 \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0))) \times (R(-\theta)\mathbf{B}(\mathbf{x}^0)). \end{aligned} \quad (26) $$
We pursue now with the computation of the second term, $T^2$, in the formula of $\tilde{\alpha}^1$. Thus, by deriving in $(\mathbf{x}, \mathbf{v})$ in (23) we obtain

$$ \nabla_{(\mathbf{x},\mathbf{v})}(\theta \tilde{\mathbf{A}}^0)(\theta, \mathbf{x}^0, \mathbf{v}^0) = \begin{pmatrix} O_3 & \mathcal{R}(\theta) \\ -\mathcal{R}(-\theta)\nabla\mathbf{E}(\mathbf{x}^0) + \nabla_{\mathbf{x}}(\mathbf{v}^0 \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0))) & M_{-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0)} \end{pmatrix} $$

and then, by (22), we compute $T^2 := -\nabla_{(\mathbf{x},\mathbf{v})}(\theta \tilde{\mathbf{A}}^0)(\theta, \mathbf{x}^0, \mathbf{v}^0)\{\tilde{\mathbf{a}}^0(\mathbf{x}^0, \mathbf{v}^0)\}$ as follows

$$ \begin{aligned} T_{\mathbf{x}}^2 &= -\mathcal{R}(\theta)\mathbf{E}(\mathbf{x}^0)_{\parallel} - \mathcal{R}(\theta)(\mathbf{v}^0 \times \mathbf{B}(\mathbf{x}^0)_{\parallel}) = -\mathcal{R}(\theta)(\mathbf{v}^0 \times \mathbf{B}(\mathbf{x}^0)_{\parallel}), \\ T_{\mathbf{v}}^2 &= \mathcal{R}(-\theta)\nabla\mathbf{E}(\mathbf{x}^0)\mathbf{v}^0_{\parallel} - \nabla_{\mathbf{x}}(\mathbf{v}^0 \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0)))\mathbf{v}^0_{\parallel} \\ &\quad - (\mathbf{E}(\mathbf{x}^0)_{\parallel} + \mathbf{v}^0 \times \mathbf{B}(\mathbf{x}^0)_{\parallel}) \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^0)). \end{aligned} \quad (27) $$

Hence, $\tilde{\alpha}^1$ is obtained as $\tilde{\alpha}^1 = T^1 + T^2$. Before concluding with this formula we note that both $T_{\mathbf{v}}^1$ and $T_{\mathbf{v}}^2$ have a term of the type $\nabla_{\mathbf{x}}(\mathbf{v} \times \mathbf{B}(\mathbf{x}))\mathbf{u}$. Thus, using the fact that for any vectors $\mathbf{u}$ and $\mathbf{v}$ and any field $\mathbf{B} \equiv \mathbf{B}(\mathbf{x})$ we have

$$ \nabla_{\mathbf{x}}(\mathbf{v} \times \mathbf{B}(\mathbf{x}))\mathbf{u} = \mathbf{v} \times (\mathbf{u} \cdot \nabla_{\mathbf{x}})\mathbf{B}(\mathbf{x}), \quad \text{where } (\mathbf{u} \cdot \nabla_{\mathbf{x}})\mathbf{B} := \sum_{i} \mathbf{u}_i \frac{\partial \mathbf{B}}{\partial \mathbf{x}_i}, $$

we can write the final form of $\tilde{\alpha}^1$ as follows:

$$ \tilde{\alpha}_{\mathbf{x}}^{1}(\theta, \mathbf{x}^{0}, \mathbf{v}^{0}, \mathbf{x}^{1}, \mathbf{v}^{1}) = R(\theta)\mathbf{v}^{1} + \mathcal{R}(\theta)\mathbf{E}(\mathbf{x}^{0}) + R(\theta)\mathbf{v}^{0} \times \mathcal{R}(\theta)\mathbf{B}(\mathbf{x}^{0}) - \mathcal{R}(\theta)(\mathbf{v}^{0} \times \mathbf{B}(\mathbf{x}^{0}))_{\parallel}, \quad (28) $$

$$ \begin{aligned} \tilde{\alpha}_{\mathbf{v}}^{1}(\theta, \mathbf{x}^{0}, \mathbf{v}^{0}, \mathbf{x}^{1}, \mathbf{v}^{1}) = & R(-\theta)\nabla\mathbf{E}(\mathbf{x}^{0})[\mathbf{x}^{1} + \mathcal{R}(\theta)\mathbf{v}^{0}] + \mathcal{R}(-\theta)\nabla\mathbf{E}(\mathbf{x}^{0})\mathbf{v}^{0}_{\parallel} \\ & + (\mathbf{v}^{1} - \mathcal{R}(-\theta)\mathbf{E}(\mathbf{x}^{0}) + \mathbf{v}^{0} \times (-\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^{0}))) \times (R(-\theta)\mathbf{B}(\mathbf{x}^{0})) \\ & + (\mathbf{E}(\mathbf{x}^{0})_{\parallel} + \mathbf{v}^{0} \times \mathbf{B}(\mathbf{x}^{0})_{\parallel}) \times \mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^{0}) \\ & + \mathbf{v}^{0} \times R(-\theta)[((\mathbf{x}^{1} + \mathcal{R}(\theta)\mathbf{v}^{0}) \cdot \nabla)\mathbf{B}(\mathbf{x}^{0})] \\ & + \mathbf{v}^{0} \times (\mathbf{v}^{0}_{\parallel} \cdot \nabla)(\mathcal{R}(-\theta)\mathbf{B}(\mathbf{x}^{0})). \end{aligned} \quad (29) $$

Finally, the last step for finding equations (19) and (20) is achieved by integrating (28)-(29) in $\theta$ over $[0, 2\pi]$. Indeed, if $(\mathbf{y}^0(t), \mathbf{u}^0(t))$ is the solution to problem (18), then $(\mathbf{y}^1(t), \mathbf{u}^1(t))$ is solution to the problem

$$ \left\{ \begin{aligned} \frac{d\mathbf{y}^1}{dt} &= \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}_x^1(\theta, \mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)d\theta, & \mathbf{y}^1(s) &= 0, \\ \frac{d\mathbf{u}^1}{dt} &= \frac{1}{2\pi} \int_0^{2\pi} \tilde{\alpha}_v^1(\theta, \mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)d\theta, & \mathbf{u}^1(s) &= 0. \end{aligned} \right. $$

These integrations can be simply computed with a symbolic software, e.g., with Maple. By an usual implementation of the operations with the matrices $R(\theta)$ and $\mathcal{R}(\theta)$ and of the cross product we obtain, for example, that

$$ \frac{1}{2\pi} \int_{0}^{2\pi} R(\theta) \mathbf{u}^{0} \times \mathcal{R}(\theta) \mathbf{B}(\mathbf{y}^{0}) - \mathcal{R}(\theta) (\mathbf{u}^{0} \times \mathbf{B}(\mathbf{y}^{0}))_{\parallel} \, d\theta = \begin{pmatrix} \mathbf{u}_{y}^{0} \mathbf{B}_{y}(\mathbf{y}^{0}) + \mathbf{u}_{z}^{0} \mathbf{B}_{z}(\mathbf{y}^{0}) \\ \mathbf{u}_{y}^{0} \mathbf{B}_{x}(\mathbf{y}^{0}) + \mathbf{u}_{x}^{0} \mathbf{B}_{y}(\mathbf{y}^{0}) \\ \mathbf{u}_{z}^{0} \mathbf{B}_{x}(\mathbf{y}^{0}) + \mathbf{u}_{x}^{0} \mathbf{B}_{z}(\mathbf{y}^{0}) \end{pmatrix} $$
and thus, integrating (28) implies immediately the equation for $\mathbf{y}^1$ in (19). The same approach for integrating the terms in (29) leads to the equation defined in (19)-(20). □

**Remark 2.2.** Using the first two terms in expansion (15) and the results of Section 2.1, we obtain for small $\epsilon$ the following approximations of the solution $(\mathbf{x}_{\epsilon}(t), \mathbf{v}_{\epsilon}(t))$ to problem (1):

(i) at the zero-th order

$$G_{\epsilon}^{0}(t) = \begin{pmatrix} \mathbf{y}^{0}(t) \\ R(\frac{t-s}{\epsilon})\mathbf{u}^{0}(t) \end{pmatrix}, \qquad (30)$$

where $(\mathbf{y}^0, \mathbf{u}^0)$ are given by (18),

(ii) at the first order

$$G_{\epsilon}^{1}(t) = \begin{pmatrix} \mathbf{y}^{0}(t) \\ R(\frac{t-s}{\epsilon})\mathbf{u}^{0}(t) \end{pmatrix} + \epsilon \left\{ \begin{pmatrix} \mathbf{y}^{1}(t) \\ R(\frac{t-s}{\epsilon})\mathbf{u}^{1}(t) \end{pmatrix} + \begin{pmatrix} \mathcal{R}(\frac{t-s}{\epsilon})\mathbf{u}^{0}(t) \\ \mathcal{R}(\frac{t-s}{\epsilon})\mathbf{E}(\mathbf{y}^{0}(t)) + R(\frac{t-s}{\epsilon})\mathbf{u}^{0}(t) \times \mathcal{R}(\frac{t-s}{\epsilon})\mathbf{B}(\mathbf{y}^{0}(t)) \end{pmatrix} \right\}, \quad (31)$$

where $(\mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)$ are given by (18)-(19)-(20) and the rotation matrices $R(\cdot)$ and $\mathcal{R}(\cdot)$ are defined in (5).

**Remark 2.3.** Unlike the zero-order model, the first-order approximation captures

* the additional force due to the term $\mathbf{B}$, in the axial direction, that was alluded to in Introduction. This force is visible at the slow scale of the parallel bounce.

* the electric drift $\mathbf{E} \times \mathbf{B}_0$ term in the equation for $\mathbf{y}^1$ in (19) which can be seen at the slowest time scale.

**Remark 2.4.** The first-order system (18)-(19)-(20) is source-free, meaning that in the enlarged variable $Y = (\mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)$, the system writes under the form

$$\frac{dY}{dt} = F(Y) \quad \text{where } F : \mathbb{R}^{12} \to \mathbb{R}^{12} \text{ satisfies } \nabla_Y \cdot F = 0. \qquad (32)$$

Indeed, denoting the components of $F$ by $F_{\mathbf{y}^0}, F_{\mathbf{u}^0}, F_{\mathbf{y}^1}$, and $F_{\mathbf{u}^1}$, it is easy to see that we have

$$\nabla_{\mathbf{y}^0} \cdot F_{\mathbf{y}^0} = \nabla_{\mathbf{u}^0} \cdot F_{\mathbf{u}^0} = \nabla_{\mathbf{y}^1} \cdot F_{\mathbf{y}^1} = \nabla_{\mathbf{u}^1} \cdot F_{\mathbf{u}^1} = 0.$$

Thus, by Liouville's theorem, the system (18)-(19)-(20) conserves the volumes in the enlarged phase space $(\mathbf{y}^0, \mathbf{u}^0, \mathbf{y}^1, \mathbf{u}^1)$.
# 3 The numerical approach

In this section we describe the time-stepping methods for solving the initial model in (1), the zero-order approximation model in (18), and the first-order model in (18)-(19)-(20). Since the three systems are source-free, we use numerical methods which conserve the volumes in the corresponding spaces, the so-called volume-preserving integrators (see [4]). Following the approach in [1] (see also [6, 7]), we can derive simple and efficient volume-preserving schemes of high order by, first, splitting the problem into several subproblems and second, by compositions of the subproblems’ solvers. The advantage of the splitting is that the subproblems can be solved analytically, unlike the initial problem. Additionally, the specific form of the composition of the subsolvers ensures the time-symmetry property, satisfied by the solution of the problem (1).

## 3.1 The volume-preserving schemes

First, we use standard volume-preserving integrators for solving (1) and (18), then we propose a volume-preserving method for solving (18)-(19)-(20). To this end, we introduce the following notation. If for any positive integer $n$, the vector $y_s \in \mathbb{R}^n$ and the vector field $F: \mathbb{R}^n \to \mathbb{R}^n$ are given, we consider the generic initial-value problem

$$ \frac{dY}{dt} = F(Y), \quad Y(s) = y_s. \qquad (33) $$

Then, for any $t > 0$ we denote by $\varphi(t) := Y(s + t)$, the solution at time $s + t$ starting from the initial value $y_s$ at time $s$. In this setting, we classically decompose the vector field of equation (1) as $F = F_1 + F_2 + F_3$, where

$$ F_1(\mathbf{x}, \mathbf{v}) = \begin{pmatrix} \mathbf{v} \\ \mathbb{O}_3 \end{pmatrix}, F_2(\mathbf{x}, \mathbf{v}) = \begin{pmatrix} \mathbb{O}_3 \\ \mathbf{E}(\mathbf{x}) \end{pmatrix}, \text{ and } F_3(\mathbf{x}, \mathbf{v}) = \begin{pmatrix} \mathbb{O}_3 \\ \mathbf{v} \times \left( \frac{1}{\varepsilon} \mathbf{B}_0 + \mathbf{B}(\mathbf{x}) \right) \end{pmatrix}. $$

The solutions, denoted by $\varphi_1, \varphi_2, \varphi_3$, of the systems (33) corresponding to $F_1, F_2, F_3$ above are easy to derive ([6, 7]). Then, if $h$ is a time step, we consider the scheme defined by

$$ G_h^4 := G_{\alpha h}^2 \circ G_{\beta h}^2 \circ G_{\alpha h}^2, \text{ where } G_h^2 := \varphi_1(h/2) \circ \varphi_2(h/2) \circ \varphi_3(h) \circ \varphi_2(h/2) \circ \varphi_1(h/2), \quad (34) $$

and where $\alpha = 1/(2 - 2^{1/3}), \beta = 1 - 2\alpha$. This scheme is volume-preserving, time-symmetric and approximates the solution of (1) to order 4 ([1]). Similarly, for the limit model in (18), we decompose the total field as $F = F_1 + F_2$ with

$$ F_1(\mathbf{y}^0, \mathbf{u}^0) = \begin{pmatrix} \mathbf{u}_x^0 \\ \mathbb{O}_2 \\ \mathbf{E}_x(\mathbf{y}^0) \\ \mathbb{O}_2 \end{pmatrix} \quad \text{and} \quad F_2(\mathbf{y}^0, \mathbf{u}^0) = \begin{pmatrix} \mathbb{O}_4 \\ \mathbf{B}_x(\mathbf{y}^0)\mathbf{u}_z^0 \\ -\mathbf{B}_x(\mathbf{y}^0)\mathbf{u}_y^0 \end{pmatrix}. $$

Denoting by $\varphi_1$ and $\varphi_2$ the solutions of the associated systems (that will be derived in Section 3.2) when the field $\mathbf{E}$ is given by (2)), we can use as before, the four-order time-stepping scheme

$$ G_h^4 := G_{\alpha h}^2 \circ G_{\beta h}^2 \circ G_{\alpha h}^2, \text{ where } G_h^2 := \varphi_1(h/2) \circ \varphi_2(h) \circ \varphi_1(h/2). \qquad (35) $$
Next, we develop the numerical method for solving the first-order model (18)-(19)-(20). This is done by means of a splitting scheme which preserves the volume in the enlarged space introduced in Remark 2.4. However, the issue of this approach is that the splitting can be performed in various ways. This difficulty is accentuated when the total force field is composed of a big number of terms.

Our choice of splitting is based on two facts: first, to be able to solve analytically the subproblems and second, to mimic the structural difficulty of the original model; more precisely, if $\mathbf{B} \equiv 0$ in equation (1) (the case of an ideal Penning trap), the motion in the $\vec{e}_x$ direction is decoupled from the motion in the perpendicular plane to $\vec{e}_x$ (see [8]). Even if this behaviour is not true when $\mathbf{B} \neq 0$, we choose this property to split, as detailed below. We thus consider the splitting method where the vector field $F$ introduced in (32) is splitted into four components. More precisely, we write $F = \sum_{i=1}^4 F_i$, where the functions $F_i : \mathbb{R}^{12} \to \mathbb{R}^{12}$ are defined by

$$F_1(Y) := \begin{pmatrix} \mathbf{u}_x^0 \\ \mathcal{O}_2 \\ \mathbf{E}_x(\mathbf{y}^0) \\ \mathcal{O}_8 \end{pmatrix}, F_2(Y) := \begin{pmatrix} \mathcal{O}_6 \\ \mathbf{u}_x^1 \\ \mathcal{O}_2 \\ \nabla \mathbf{E}_x(\mathbf{y}^0) \cdot \mathbf{y}^1 \end{pmatrix} + \begin{pmatrix} \mathcal{O}_6 \\ \mathbf{u}_y^0 \mathbf{B}_y(\mathbf{y}^0) + \mathbf{u}_z^0 \mathbf{B}_z(\mathbf{y}^0) \\ \mathcal{O}_2 \\ \mathcal{J}_x(\mathbf{y}^0, \mathbf{u}^0) \\ \mathcal{O}_2 \end{pmatrix}, \quad (36)$$

$$F_3(Y) := \begin{pmatrix} \mathcal{O}_4 \\ \mathbf{B}_x(\mathbf{y}^0)\mathbf{u}_z^0 \\ -\mathbf{B}_x(\mathbf{y}^0)\mathbf{u}_y^0 \\ \mathcal{O}_6 \end{pmatrix}, \quad F_4(Y) := \begin{pmatrix} \mathcal{O}_7 \\ \mathbf{E}_z(\mathbf{y}^0) + \mathbf{u}_y^0 \mathbf{B}_x(\mathbf{y}^0) + \mathbf{u}_x^0 \mathbf{B}_y(\mathbf{y}^0) \\ -\mathbf{E}_y(\mathbf{y}^0) + \mathbf{u}_z^0 \mathbf{B}_x(\mathbf{y}^0) + \mathbf{u}_x^0 \mathbf{B}_z(\mathbf{y}^0) \\ 0 \\ \mathbf{B}_x(\mathbf{y}^0)\mathbf{u}_z^1 + \mathbf{u}_z^0 \nabla \mathbf{B}_x(\mathbf{y}^0) \cdot \mathbf{y}^1 + \mathcal{J}_y(\mathbf{y}^0, \mathbf{u}^0) \\ -\mathbf{B}_x(\mathbf{y}^0)\mathbf{u}_y^1 - \mathbf{u}_y^0 \nabla \mathbf{B}_x(\mathbf{y}^0) \cdot \mathbf{y}^1 + \mathcal{J}_z(\mathbf{y}^0, \mathbf{u}^0) \end{pmatrix}, \quad (37)$$

and where $\mathcal{J}$ is given in (20). The differential equations corresponding to the vector field $F_1 + F_2$ describe the evolution in time of $(\mathbf{y}_x^0, \mathbf{u}_x^0, \mathbf{y}_z^1, \mathbf{u}_z^1)$, while those corresponding to $F_3 + F_4$ describe the evolution of $(\mathbf{y}_\perp^0, \mathbf{u}_\perp^0, \mathbf{y}_\perp^1, \mathbf{u}_\perp^1)$. Writing the initial-value problem in (33) for each $(F_i)_{i \in \{1, \dots, 4\}}$, solutions of which are denoted by $\varphi_i(t)$, we obtain from [1] the following result.

**Proposition 3.1.** Let the vector field $F$ of the system (18)-(19)-(20) be decomposed as $F = \sum_{i=1}^4 F_i$, where $(F_i)_{i \in \{1, \dots, 4\}}$ are given in (36)-(37). Then, if $h$ denotes a generic time step, we construct the following time-stepping method

$$G_h^4 := G_{\alpha h}^2 \circ G_{\beta h}^2 \circ G_{\alpha h}^2, \text{ where} \quad (38)$$

$$G_h^2 := \varphi_1(h/2) \circ \varphi_2(h/2) \circ \varphi_3(h/2) \circ \varphi_4(h) \circ \varphi_3(h/2) \circ \varphi_2(h/2) \circ \varphi_1(h/2), \quad (39)$$

and where $\alpha = 1/(2 - 2^{1/3})$ and $\beta = 1 - 2\alpha$. This scheme is volume-preserving, time-symmetric and approximates the solution of (18)-(19)-(20) to order 4.

**Remark 3.2.** Following [1], by compositions of low-order schemes we can build a time-stepping method of arbitrarily even high-order. However, for simplicity, in this work we restricted to four-order schemes.

In the next section we give the analytic solutions of the initial-value problems associated to $(F_i)_{i \in (36)-(37)}$, in the case of the electric field $\mathbf{E}$ in (2).
## 3.2 The analytic solutions of the first-order subsystems

In this section we make further explicit the model in (18)-(19)-(20), by considering the electric field given in (2), whereas the magnetic field **B** is kept general. In this case, we have that the jacobian matrix of **E** and, for example, the electric term in (20) simplify to

$$ \nabla \mathbf{E}(\mathbf{x}) = \begin{pmatrix} -c & 0 & 0 \\ 0 & c/2 & 0 \\ 0 & 0 & c/2 \end{pmatrix} \quad \text{and} \quad \mathcal{J}_{\delta}(\mathbf{y}^0, \mathbf{u}^0) = \frac{c}{2} \begin{pmatrix} 0 \\ -\mathbf{u}_z^0 \\ \mathbf{u}_y^0 \end{pmatrix}, $$

respectively. In the sequel, we specify the initial-value problems corresponding to the fields $(F_i)_{i=1,...,4}$ in (36)-(37). First, by neglecting the equations of the components which are constant in time, the equations for $F_1$ and $F_2$ become

$$ \frac{d}{dt} \begin{pmatrix} \mathbf{y}_x^0 \\ \mathbf{u}_x^0 \end{pmatrix} = \begin{pmatrix} \mathbf{u}_x^0 \\ -c\mathbf{y}_x^0 \end{pmatrix} \quad \text{and} \quad \frac{d}{dt} \begin{pmatrix} \mathbf{y}_x^1 \\ \mathbf{u}_x^1 \end{pmatrix} = \begin{pmatrix} \mathbf{u}_x^1 \\ -c\mathbf{y}_x^1 \end{pmatrix} + \begin{pmatrix} k_y \\ k_u \end{pmatrix}, \quad (40) $$

where

$$ k_y := \mathbf{u}_y^0 \mathbf{B}_y(\mathbf{y}^0) + \mathbf{u}_z^0 \mathbf{B}_z(\mathbf{y}^0) \quad \text{and} \quad k_u := \mathcal{J}_x(\mathbf{y}^0, \mathbf{u}^0) $$

are constants in time thanks to the splitting approach. Second, by the same procedure, the equations for $F_3$ and $F_4$ become

$$ \frac{d}{dt} \begin{pmatrix} \mathbf{u}_y^0 \\ \mathbf{u}_z^0 \end{pmatrix} = \begin{pmatrix} r\mathbf{u}_y^0 \\ -r\mathbf{u}_z^0 \end{pmatrix} \quad \text{and} \quad \frac{d}{dt} \begin{pmatrix} \mathbf{y}_y^1 \\ \mathbf{y}_z^1 \\ \mathbf{u}_y^1 \\ \mathbf{u}_z^1 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ r\mathbf{u}_z^1 + \mathbf{y}^1 \cdot \mathbf{\Lambda}_y \\ -r\mathbf{u}_y^1 - \mathbf{y}^1 \cdot \mathbf{\Lambda}_z \end{pmatrix} + \begin{pmatrix} \ell_y \\ \ell_z \\ m_y \\ m_z \end{pmatrix}, \quad (41) $$

where

$$ r := \mathbf{B}_x(\mathbf{y}^0), \quad \mathbf{\Lambda}_y := \mathbf{u}_z^0 \nabla \mathbf{B}_x(\mathbf{y}^0), \quad \mathbf{\Lambda}_z := \mathbf{u}_y^0 \nabla \mathbf{B}_x(\mathbf{y}^0), $$

$$ \ell_y := \mathbf{E}_z(\mathbf{y}^0) + \mathbf{u}_y^0 \mathbf{B}_x(\mathbf{y}^0) + \mathbf{u}_x^0 \mathbf{B}_y(\mathbf{y}^0), \quad \ell_z := -\mathbf{E}_y(\mathbf{y}^0) + \mathbf{u}_z^0 \mathbf{B}_x(\mathbf{y}^0) + \mathbf{u}_x^0 \mathbf{B}_z(\mathbf{y}^0), $$

$$ m_y := \mathcal{J}_y(\mathbf{y}^0, \mathbf{u}^0) \quad \text{and} \quad m_z := \mathcal{J}_z(\mathbf{y}^0, \mathbf{u}^0) $$

are constants in time. Equations (40)-(41) are easy to solve, their analytic solutions are given below. For any time $t$ and any time step $h > 0$ we have the solutions of the systems in (40):

$$ \begin{pmatrix} \mathbf{y}_x^0(t) \\ \mathbf{u}_x^0(t) \end{pmatrix} = \begin{pmatrix} \cos(\sqrt{ch}) & \sin(\sqrt{ch})/\sqrt{c} \\ -\sqrt{c}\sin(\sqrt{ch}) & \cos(\sqrt{ch}) \end{pmatrix} \begin{pmatrix} \mathbf{y}_x^0(t-h) \\ \mathbf{u}_x^0(t-h) \end{pmatrix} \quad \text{and} $$

$$ \begin{pmatrix} \mathbf{y}_x^1(t) \\ \mathbf{u}_x^1(t) \end{pmatrix} = \begin{pmatrix} \cos(\sqrt{ch}) & \sin(\sqrt{ch})/\sqrt{c} \\ -\sqrt{c}\sin(\sqrt{ch}) & \cos(\sqrt{ch}) \end{pmatrix} \begin{pmatrix} \mathbf{y}_x^1(t-h) - k_u/c \\ \mathbf{u}_x^1(t-h) + k_y \end{pmatrix} + \begin{pmatrix} k_u/c \\ -k_y \end{pmatrix}, $$

respectively. As for the systems in (41) we can compute

$$ \begin{pmatrix} \mathbf{u}_y^0(t) \\ \mathbf{u}_z^0(t) \end{pmatrix} = \begin{pmatrix} \cos(rh) & \sin(rh) \\ -\sin(rh) & \cos(rh) \end{pmatrix} \begin{pmatrix} \mathbf{u}_y^0(t-h) \\ \mathbf{u}_z^0(t-h) \end{pmatrix} \quad \text{and} $$

$$ \begin{pmatrix} \mathbf{y}_y^1(t) \\ \mathbf{y}_z^1(t) \end{pmatrix} = h \begin{pmatrix} \ell_y \\ \ell_z \end{pmatrix} + \begin{pmatrix} \mathbf{y}_y^1(t-h) \\ \mathbf{y}_z^1(t-h) \end{pmatrix}, $$
while for deriving the formulas of $\mathbf{u}_y^1$ and $\mathbf{u}_z^1$, we take into account that $\mathbf{y}_x^1$ is constant and $\mathbf{y}_y^1$ and $\mathbf{y}_z^1$ are given above. Thus, for some constants in time $n_y, n_z, p_y, p_z$, we can rewrite the differential equations in (41) under the form

$$ \frac{d}{dt} \begin{pmatrix} \mathbf{u}_y^1 \\ \mathbf{u}_z^1 \end{pmatrix} = r \begin{pmatrix} \mathbf{u}_z^1 \\ -\mathbf{u}_y^1 \end{pmatrix} + h \begin{pmatrix} n_y \\ n_z \end{pmatrix} + \begin{pmatrix} p_y \\ p_z \end{pmatrix} $$

which solution is simple to compute.

Then, we implement these integrators, in the cases where the magnetic part $\mathbf{B}$ is given by (3) and by (4). In particular, all the terms in (18)-(19)-(20) are explicit and the equations are much simplified.

## 4 Numerical results. Application to Penning traps

In this section we detail the implementation of the numerical methods for the three models, in order to illustrate the approximation properties previously described. First, we recall that the starting test case is that of an ideal Penning trap [8], a device for storing charged particles by means of an homogeneous magnetic field (the $\mathbf{B}_0/\varepsilon$ term in (1)) and a static and inhomogeneous quadrupole electric field, given by (2). In such a trap, a particle performs three motions with characteristic frequencies: a modified cyclotron motion, at frequency around $1/\varepsilon$, an axial motion, at frequency $\sqrt{c}$, and a magnetron motion or the $\mathbf{E} \times \mathbf{B}_0$ drift, at frequency around $\varepsilon$.

We emphasize that in practice the magnetic field involves perturbations, described by the term $\mathbf{B}$ in (1). The test cases of the magnetic field bottle and the asymmetric magnetic field are classic for studying the problem of a realistic Penning trap. Below, we perform numerical simulations by using the expressions in (3) and (4) of the perturbation term $\mathbf{B}$, as in [7]. We fixed the values $c = 2$ for the electric field in (2) and $k = 4$ for the amplitude of the magnetic perturbation. Clearly, increasing $k$ takes the particle further away from the trajectory in an ideal Penning trap (see [7, Fig. 3]). Then, increasing the parameter $c$ makes both, the axial time period and the (large) magnetron period, smaller. Therefore, when $c$ grows, the particle travels a greater distance in a fixed time interval.

In the following we consider the initial condition $\mathbf{x} = (1, 0, 1)^T, \mathbf{v} = (1, 1, 1)^T$. For the time step of the reduced models, referred to in the sequel as the coarse time step, we need to solve the period $T_a \approx 2\pi/\sqrt{c}$ of the axial slow motion of the particle (see Fig. 3 for an illustration). We took 9 points for solving this period which results in a coarse time step of $\Delta t \approx 0.5$. As for the stiff model (1), we have to accurately resolve the smallest time scale and therefore, we took the fine time step $\delta t = 2\pi\varepsilon/40$. However, when assessing the efficiency of the numerical methods we also consider fine solvers using 10 or 20 points for solving the cyclotron period. We recall that for the three models we use volume-preserving integrators of order 4. From now on, we call (i) the numerical scheme that discretize equation (1) the fine integrator, (ii) the scheme for the first-order system in (18)-(19)-(20) the $L_1$ integrator, and (iii) the scheme for the zero-order model in (18) the $L_0$ integrator.
## 4.1 Numerical accuracy and efficiency of the first-order integrator

In the sequel, we define the relative error at time $t_n = n\Delta t$ with respect to the fine solver $\mathcal{F}$ of a reduced model integrator $\mathcal{G}$ by

$$ \text{Err}(t_n) = \frac{\|\mathcal{G}(t_n) - \mathcal{F}(t_n)\|_1}{\|\mathcal{F}(t_n)\|_1}, \qquad (42) $$

where $\| \cdot \|_1$ stands for the $\ell_1$ norm in $\mathbb{R}^6$. First, we plot in Fig. 1 the time evolution of the relative error defined in (42) of the $L_1$ integrator with respect to the fine one, until $T = 10$. We observe, as expected, the property of smaller errors with smaller $\varepsilon$. Then, we can see that in small final times, of order 1, the first-order model is a very good approximation, as the theory predicts, whereas, in larger times the quality of the approximation is not preserved for all the values of $\varepsilon$.

Next, we remark that the error of the $L_1$ integrator is much smaller than that of the $L_0$ solver (see Fig. 2). The reasons are described in Remark 2.3. To illustrate the first point of this remark, we plot in Fig. 3 the evolution in time of the position and of the velocity in the axial ($\vec{e}_x$) direction. We can see that the $L_0$ integrator does not capture the correct period of this motion, unlike the $L_1$ integrator.

Then, we also implemented an explicit Runge-Kutta scheme of order 4 for the reduced first-order model<fnref n="2" />, in order to pinpoint that such a scheme does not provide accurate results. Indeed, it is known [4] that volume-preserving methods are preferred, in order to conserve the geometric properties of the system. Typically, we first checked that the conservation in time of the energy is respected for the volume-preserving method, unlike the RK integrator. Moreover, in our cases, the quantities $\mathbf{v}_y^2 + \mathbf{v}_z^2$ and $\mathbf{x}_x^2 + \mathbf{v}_x^2$ are not invariant in time but slowly varying. Thus, we next verified that the volume-preserving integrator keeps these quantities bounded in time, whereas the RK method spuriosly vanishes them (simulations not shown here). In addition to the geometric invariants, the numerical error is known to accumulate more slowly with the volume-preserving integrator. We checked this property, in Fig. 2.

Finally, we investigate the efficiency of the reduced models in terms of computation time, only in the case of the magnetic bottle. With this aim, we first measure the simulation time of three fine solvers together with their corresponding accuracy. More precisely, we compute the relative error of a fine solver against a very fine solution, called hereafter the reference solution, since an analytic one is not available for the system (1)-(2)-(3). We define the reference solution as the one computed with time step $2\pi\varepsilon/400$. Second, we compute the run time of the $L_0$ and $L_1$ solvers. We report the outcomes of these simulations in Tab. 1. The last line of the table corresponds to the integer number of the fine time steps in one coarse time step, which is the theoretical gain of a reduced model over the stiff one. The measured time ratios reported in lines 4 and 5 are clearly smaller than the theoretical ones, since

$$ \text{Time}(\Delta t(L_1)) > \text{Time}(\Delta t(L_0)) > \text{Time}(\delta t(\text{fine solver})), $$

<footnote n="2">²In the illustrations we call this method the RK $L_1$ integrator.</footnote>
where we denoted by $Time(\delta t(\mathcal{S}))$ the computation time of one time step of a solver $\mathcal{S}$. These inequalities hold because of the much simpler form of the equation in the stiff model when compared to the reduced one. In addition, the number of stages in the numerical algorithm of the $L_1$ integrator is larger than the one in the $L_0$ integrator.

The results in Tab. 1 show that the $L_1$ integrator is 48 times faster than the fine integrator of the initial stiff model with $\varepsilon = 10^{-2}$, in the case of a target relative error of the fine solver of $10^{-3}$ for the final time $T = 80$. If an accuracy of order $10^{-6}$ is aimed for the stiff solution with $\varepsilon = 10^{-3}$, the gain factor of the first-order model is 1739. Note that the $L_0$ integrator is clearly slightly faster than the $L_1$ integrator, but it provides a much lower accuracy, as shown at the beginning of this section.

<figure><img src="image_2.png" /><figcaption>Figure 1: The relative error of the first-order model, for several values of $\varepsilon$. Left: The magnetic bottle case, Right: The asymmetric magnetic field case.</figcaption></figure>

## 4.2 Parareal long-time simulations

In the previous section, we observed that the first-order model provides a good approximation in short times, of order 10 at most. However, the accuracy of this approximation is lost in large times, of order $1/\varepsilon$. Thus, we can see in Fig. 4 that for the largest values of $\varepsilon$, the first-order model is not reliable anymore at time $T = 80$, which corresponds to a large fraction of the magnetron cycle.

To overcome this issue, we next show that within the Parareal framework [11], we can retrieve accurate approximation in large times. In the following, we briefly recall the Parareal algorithm in the framework of the initial-value problem (1) to be solved over the interval $[0, T]$. We decompose this interval into $N$ uniform time intervals $[T_n, T_{n+1}]$ for $n \in \{0, \dots, N-1\}$. Let $\mathcal{F}(T_{n+1}, T_n, U_n)$ denote the value of a fine solver providing an accurate approximation of the solution of (1) at time $T_{n+1}$ with the initial condition $U_n$ at time $T_n$. Let $\mathcal{G}(T_{n+1}, T_n, U_n)$ denote the value of a coarse solver providing a less accurate approximation to the same solution. In our framework, the computational cost of the $\mathcal{F}$ term is relatively high since its time step is constrained by the small parameter $\varepsilon$, while the $\mathcal{G}$ solver should be chosen such that its cost is much lower than the one of the fine solver. Following the strategy introduced in [3], which goes back to [12], we use the first-order reduced model for the coarse solving. Thus, the Parareal algorithm consists in computing a
<figure><img src="image_3.png" /><figcaption>Figure 2: Volume-preserving vs. Runge-Kutta 4 methods and first-order vs. zero-order models, for two values of $\epsilon$. Top: The magnetic bottle case, Bottom: The asymmetric magnetic field case.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: The time evolution of the position and velocity in the axial direction in the magnetic bottle case, for $\epsilon = 0.01$.</figcaption></figure>

sequence $(U_n^k)_{n,k}$ of approximations of $(\mathbf{x}_\epsilon, \mathbf{v}_\epsilon)(T_n)$ for every $k \in \{0, 1, \dots\}$ in the following way: the initial approximation $(U_n^0)_n$ is computed sequentially using the coarse solver

$$ \forall n \in \{0, \dots, N-1\} \quad U_{n+1}^0 = \mathcal{G}(T_{n+1}, T_n, U_n^0), \quad \text{with } U_0^0 = (\mathbf{x}, \mathbf{v}) \qquad (43) $$
<table><tr><td>Nb. of points in P = 2πε</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Error(fine solver)</td><td>4.2543 · 10<sup>-3</sup></td><td>2.816 · 10<sup>-4</sup></td><td>1.77 · 10<sup>-5</sup></td></tr><tr><td>Time(fine solver) in sec.</td><td>1.26 · 10<sup>-2</sup></td><td>2.37 · 10<sup>-2</sup></td><td>4.79 · 10<sup>-2</sup></td></tr><tr><td>Time(fine solver)/Time(L<sub>0</sub> solver)</td><td>60.0</td><td>122.0</td><td>234.7</td></tr><tr><td>Time(fine solver)/Time(L<sub>1</sub> solver)</td><td>47.9</td><td>93.6</td><td>178.7</td></tr><tr><td>Upper bound of Time(fine solver)/Time(reduced solver)</td><td>79</td><td>157</td><td>314</td></tr></table>

<table><caption>Table 1: The computation time ratios between the reduced and the fine integrators vs. the error of the fine integrator at the final time $T = 80$. The time step of the reduced model is $\Delta t \approx 0.5$. The case of the magnetic bottle with $\epsilon = 0.01$ (top panel) and with $\epsilon = 0.001$ (bottom panel).</caption><tr><td>Nb. of points in P = 2πε</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Error(fine solver)</td><td>4.5956 · 10<sup>-4</sup></td><td>2.989 · 10<sup>-5</sup></td><td>1.89 · 10<sup>-6</sup></td></tr><tr><td>Time(fine solver) in sec.</td><td>1.20 · 10<sup>-1</sup></td><td>2.35 · 10<sup>-1</sup></td><td>4.71 · 10<sup>-1</sup></td></tr><tr><td>Time(fine solver)/Time(L<sub>0</sub> solver)</td><td>546.1</td><td>1203.1</td><td>2310.3</td></tr><tr><td>Time(fine solver)/Time(L<sub>1</sub> solver)</td><td>486.2</td><td>916.4</td><td>1739.2</td></tr><tr><td>Upper bound of Time(fine solver)/Time(reduced solver)</td><td>786</td><td>1571</td><td>3143</td></tr></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: The relative error of the first-order model, for several values of $\epsilon$, in long time. Left: The magnetic bottle case, Right: The asymmetric magnetic field case.</figcaption></figure>

and then, for all $k = 0, 1, \dots$, the sequence $(U_n^k)_{n,k}$ is updated by

$$U_{n+1}^{k+1} = \mathcal{G}(T_{n+1}, T_n, U_n^{k+1}) + \mathcal{F}(T_{n+1}, T_n, U_n^k) - \mathcal{G}(T_{n+1}, T_n, U_n^k) \quad \forall n \in \{0, \dots, N-1\}, \quad (44)$$

starting from $U_0^{k+1} = (\mathbf{x}, \mathbf{v})$. In this iteration, the terms $\mathcal{F}(T_{n+1}, T_n, U_n^k)$ have the most significant computational cost, and therefore they have to be computed in parallel over the number $N$ of time windows. In this setting, the Parareal efficiency mainly depends on the number of iterations, which should be significantly smaller than $N$, and on the cost ratio between the coarse and the fine solvers [13].
In the following, we recall that we denote the coarse time step by $\Delta t = T/N$ and we compute the relative error between the Parareal numerical approximation $U_n^k$ obtained at the k-th iteration and the reference numerical solution $S(t_n)$, in $L^\infty([0, T])$ as a function of the number k of Parareal iterations

$$ \text{Error}(k) = \max_{n \in \{1, \dots, N\}} \frac{\|U_n^k - S(t_n)\|_1}{\|S(t_n)\|_1}. \quad (45) $$

The criterion for stopping the Parareal iterations is when the error of the algorithm is smaller than that of the fine solver.

The first numerical experiments tackle the use of the $L_0$ integrator for the coarse propagator in the Parareal method. We performed simulations in the case $\varepsilon = 0.01$ for both expressions of magnetic field, over four time intervals with $T \in \{10, 20, 40, 80\}$ and with a fixed coarse time step $\Delta t \approx 0.5$. Thus, the number of time windows varies in $N \in \{20, 41, 81, 162\}$. The simulations (not shown here) demonstrate in both test cases and for each value of $N$ that the relative error of the Parareal algorithm is large for almost all the $N$ iterations and becomes small enough at the iteration $N-1$ or $N$. Therefore, this Parareal strategy does not allow to reduce the simulation time and makes the use of the limit model for the coarse solving ineffective. These results are justified by the inaccuracy of the limit model. Thus, we pursue the approach by using the first-order model for the coarse solver. The first numerical illustration concerns the behaviour of the Parareal errors in (45) in short-time simulations, when the parameter $\varepsilon$ decreases. We consider the fine solver with 20 points for the smallest time scale, which is a trade-off between computation time and providing small enough error (see also the speedup analysis below). The final time is $T = 1000\varepsilon$ corresponding to approximately 159 cyclotron periods for every $\varepsilon$. We display in Fig. 5 the evolution of these errors in log-scale, consistently decreasing with decreasing $\varepsilon$.

We now show the outcome of the more challenging long time simulations. We fix $\varepsilon = 0.01$ for which the first-order model gave the worst approximation in long time. As before, we keep the coarse time step fixed to $\Delta t \approx 0.5$ while the final time is gradually increased, $T \in \{10, 20, 40, 80\}$. The error in the log-scale of the Parareal method is displayed in Fig. 6 for both test cases. For example, in the simulations until $T = 80$ the target error is attained in only 8 Parareal iterations, when the number of time windows is $N = 162$.

Finally, we analyse the efficiency of the Parareal algorithm in terms of ideal speedup. To this end, we denote by $T_c$ the measured run time of the coarse solver over the total time interval $[0, T]$ and by $T_f$ the same run time for the fine solver. Then, the formula of the ideal speedup corresponding to perform the computation of the fine solver in parallel over $N$ processors is known [13]:

$$ S(N) = \frac{1}{(1 + K)\frac{T_c}{T_f} + \frac{K}{N}} \quad (46) $$

where $K$ is the number of iterations after which the error of the Parareal algorithm with respect to the reference solution is smaller than that of the fine solver. We thus report in Tables 2 and 3 the speedup of the algorithm together with the number of iterations required to achieve the target accuracy, for three fine solvers. First, we point out that, for a fixed
precision, the error of the fine solver is larger in the case of the magnetic perturbation given by (4) than that obtained with (3), consistently with the work in [7]. In addition, we emphasize that the accuracy of the fine solver given by 10 points for solving the cyclotron period is quite low for the second case. Then, we remark in both cases of magnetic field the very good speedup numbers when $\epsilon = 0.001$, which are due to the high accuracy of the first-order model that allows the Parareal method to converge in a few iterations. In the case of $\epsilon = 0.01$ the reduced model is not accurate at time $T = 80$ (recall Fig. 4) and therefore, more iterations are needed for convergence.

# 5 Conclusion

The main goal of this work was to derive a first-order reduced model to approach the solution of a stiff Newton-Lorentz equation modeling the dynamics of a particle in a realistic Penning trap, where an inhomogeneous magnetic perturbation is added to a strong constant homogeneous magnetic field. For the numerical solution, we implemented a volume-preserving method with a specific splitting technique. Then, an accuracy and efficiency analysis of the reduced models was conducted in comparison with the original model. The first-order model turned out to be effective in approximating the stiff problem, whereas the zero-th-order model did not bring a substantial information. Next, we implemented a specific version of the Parareal algorithm in order to perform accurate and robust simulations in long time. Finally, we analysed the numerical convergence and the speedup of the Parareal method in the context of two realistic magnetic perturbations.

Despite its complicated form, the first-order model has the potential to show how perturbation terms are interacting in long scales without the necessity of solving the smallest scale. As an alternative to a long numerical procedure, a theoretical investigation about the rate of convergence and the time validity of such an approximation would be interesting to undertake. In addition, such an analysis could be integrated in the Parareal approach in order to understand the efficient convergence of the Parareal algorithm over any time interval. Indeed, taking into account the favorable cost ratios between the first-order and the fine solvers, the Parareal algorithm turns out to be an effective method to perform very long simulations of stiff equations as those addressed in this work.

**Acknowledgement:** The work of the author was supported by ANR MUFFIN (ANR-19-CE46-0004).

# References

[1] K. Feng and Z-J. Shang. Volume-preserving algorithms for source-free dynamical systems. *Numerische Mathematik*, 71:451–463, 1995.

[2] E. Frénod. Application of the averaging method to the gyrokinetic plasma. *Asymptotic Analysis*, 46:1–28, 2006.
<table><tr><td>Nb. of points in P = 2πε</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Error(fine solver)</td><td>4.2543 · 10<sup>-3</sup></td><td>2.816 · 10<sup>-4</sup></td><td>1.77 · 10<sup>-5</sup></td></tr><tr><td>Nb. of Parareal iterations</td><td>7</td><td>8</td><td>8</td></tr><tr><td>Speedup</td><td>4.8</td><td>6.7</td><td>10.0</td></tr></table>

<table><caption>Table 2: Speedup of the Parareal algorithm using the $L_1$ integrator with $\Delta t \approx 0.5$ for the coarse solving, for three fine solvers, at the final time $T = 80$. The case of the magnetic bottle with $\epsilon = 0.01$ at the top and with $\epsilon = 0.001$ at the bottom.</caption><tr><td>Nb. of points in P = 2πε</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Error(fine solver)</td><td>4.5956 · 10<sup>-4</sup></td><td>2.989 · 10<sup>-5</sup></td><td>1.89 · 10<sup>-6</sup></td></tr><tr><td>Nb. of Parareal iterations</td><td>2</td><td>3</td><td>4</td></tr><tr><td>Speedup</td><td>47.7</td><td>43.4</td><td>36.0</td></tr></table>

<table><tr><td>Nb. of points in P = 2πε</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Error(fine solver)</td><td>5.4676 · 10<sup>-1</sup></td><td>3.620 · 10<sup>-2</sup></td><td>2.29 · 10<sup>-3</sup></td></tr><tr><td>Nb. of Parareal iterations</td><td>7</td><td>8</td><td>10</td></tr><tr><td>Speedup</td><td>4.7</td><td>6.7</td><td>8.1</td></tr></table>

<table><caption>Table 3: Speedup of the Parareal algorithm using the $L_1$ integrator with $\Delta t \approx 0.5$ for the coarse solving, for three fine solvers, at the final time $T = 80$. The case of the asymmetric magnetic field with $\epsilon = 0.01$ at the top and with $\epsilon = 0.001$ at the bottom.</caption><tr><td>Nb. of points in P = 2πε</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Error(fine solver)</td><td>3.900 · 10<sup>-2</sup></td><td>2.57 · 10<sup>-3</sup></td><td>1.6 · 10<sup>-4</sup></td></tr><tr><td>Nb. of Parareal iterations</td><td>1</td><td>1</td><td>2</td></tr><tr><td>Speedup</td><td>92.1</td><td>116.2</td><td>70.1</td></tr></table>
<figure><img src="image_6.png" /><figcaption>Figure 5: Convergence of the Parareal algorithm in short time simulations. The fine time step is $2\pi\epsilon/20$. Left: The magnetic bottle case, Right: The asymmetric magnetic field case.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Convergence of the Parareal algorithm in long time simulations. The fine time step is $2\pi\epsilon/20$. Left: The magnetic bottle case, Right: The asymmetric magnetic field case.</figcaption></figure>

[3] L. Grigori, S. A. Hirstoaga, V.-T. Nguyen, and J. Salomon. Reduced model-based parareal simulations of oscillatory singularly perturbed ordinary differential equations. *Journal of Computational Physics*, 436:110282, 2021.

[4] E. Hairer, C. Lubich, and G. Wanner. *Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations*. 2nd Edition, Springer, Berlin, 2006.

[5] T. Haut and B. Wingate. An asymptotic parallel-in-time method for highly oscillatory PDEs. *SIAM Journal on Scientific Computing*, 36(2):A693–A713, 2014.

[6] Y. He, Y. Sun, J. Liu, and H. Qin. Volume-preserving algorithms for charged particle dynamics. *Journal of Computational Physics*, 281:135–147, 2015.
[7] 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.

[8] M. Kretzschmar. Particle motion in a Penning trap. *European Journal of Physics*, 12(5):240, 1991.

[9] F. Legoll, T. Lelievre, and G. Samaey. A micro-macro parareal algorithm: application to singularly perturbed ordinary differential equations. *SIAM Journal on Scientific Computing*, 35(4):A1951–A1986, 2013.

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

[11] J.-L. Lions, Y. Maday, and G. Turinici. A "parareal" in time discretization of PDE's. *Comptes Rendus de l'Académie des Sciences - Series I - Mathematics*, 332:661–668, 2001.

[12] Y. Maday. Parareal in time algorithm for kinetic systems based on model reduction. *High-dimensional partial differential equations in science and engineering, CRM Proc. Lecture Notes*, 41:183–194, 2007.

[13] M. L. Minion. A hybrid parareal spectral deferred corrections method. *Communications in Applied Mathematics and Computational Science*, 5(2):265–301, 2010.

[14] A. G. Peddle, T. Haut, and B. Wingate. Parareal convergence for oscillatory pdes with finite time-scale separation. *SIAM Journal on Scientific Computing*, 41(6):A3476–A3497, 2019.

[15] C. C. Rodegheri, K. Blaum, H. Kracke, S. Kreim, A. Mooser, W. Quint, S. Ulmer, and J. Walz. An experiment for the direct determination of the *g*-factor of a single proton in a Penning trap. *New Journal of Physics*, 14(6):063011, 2012.

[16] M. Vogel. *Particle Confinement in Penning Traps*, volume 100. Springer, 2018.

[17] 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.