# Multiscale numerical schemes for the collisional Vlasov equation in the finite Larmor radius approximation regime

Anaïs Crestetto, Nicolas Crouseilles, Damien Prel

► To cite this version:

Anaïs Crestetto, Nicolas Crouseilles, Damien Prel. Multiscale numerical schemes for the collisional Vlasov equation in the finite Larmor radius approximation regime. Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal, 2023, 21, pp.1-24. 10.1137/22M1496839. hal-03690427v2

HAL Id: hal-03690427

https://hal.science/hal-03690427v2

Submitted on 13 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 Attribution 4.0 International License
# Multiscale numerical schemes for the collisional Vlasov equation in the finite Larmor radius approximation regime

Anaïs Crestetto*

Nicolas Crouseilles†

Damien Prel*

September 14, 2023

## Abstract

This work is devoted to the construction of multiscale numerical schemes efficient in the finite Larmor radius approximation of the collisional Vlasov equation. Following the paper of Bostan and Finot (2019), the system involves two different regimes, a highly oscillatory and a dissipative regimes, whose asymptotic limits do not commute. In this work, we consider a Particle-In-Cell discretization of the collisional Vlasov system which enables to deal with the multiscale characteristics equations. Different multiscale time integrators are then constructed and analysed. We prove asymptotic properties of these schemes in the highly oscillatory regime and in the collisional regime. In particular, the asymptotic preserving property towards the modified equilibrium of the averaged collision operator is recovered. Numerical experiments are then shown to illustrate the properties of the numerical schemes.

**Keywords.** Vlasov equation, Finite Larmor radius, Multiscale numerical scheme

**MSC codes.** 65M25, 65M75, 35Q83

## 1 Introduction

In this work, we consider the three dimensional collisional Vlasov model under the so-called finite Larmor radius scaling, whose unknown is the distribution function of a population of charged particles. Motivated by the magnetic confinement, in the finite Larmor radius regime, the typical perpendicular spatial length is of the same order as the Larmor radius and the parallel spatial length is much larger (see [3, 20]). When the collision effects are neglected, the finite Larmor limit leads to the so-called gyrokinetic model. However, when a collision kernel is present, the gyrokinetic model is supplemented with an averaged collision operator whose properties are different from the original one (see [4] and [2]). Indeed, under some quite general assumptions, the authors proved in [4] that the averaged (or effective) collision operator satisfies a set of modified invariants and have a modified equilibrium (compared to the non averaged version). As a consequence, the hydrodynamic limit gives rise to a new system of conservation laws called the gyrofluid model.

More precisely, we are interested in the following collisional Vlasov model satisfied by the distribution function $f(t, x, v)$

$$ \partial_t f + \frac{1}{\epsilon}(v_1 \partial_{x_1} f + v_2 \partial_{x_2} f) + v_3 \partial_{x_3} f + E \cdot \nabla_v f + \frac{1}{\epsilon}(v_2 \partial_{v_1} f - v_1 \partial_{v_2} f) = \frac{1}{\tau} Q(f), \quad (1) $$

where $(t, x, v) \in \mathbb{R}_+ \times \mathbb{R}^3 \times \mathbb{R}^3$ ($x = (x_1, x_2, x_3)$ and $v = (v_1, v_2, v_3)$), $E = (E_1, E_2, E_3) : \mathbb{R}^3 \to \mathbb{R}^3$ is an external electric field and $Q$ a collisional operator. Here the external magnetic field is supposed to be homogeneous and oriented along the $x_3$ direction. An initial condition $f(0, x, v) = f_0(x, v)$ and boundary conditions enable to complete the model.

Our main motivation is to design and analyze multiscale numerical schemes for (1), which can run for arbitrary values of $\epsilon$ (scaled cyclotronic period) and $\tau$ (Knudsen number). Let us remark that the limit $\epsilon \to 0$ ($\tau$ fixed) leads to highly-oscillatory (in time) behavior of $f$ whereas the limit $\tau \to 0$ ($\epsilon$ fixed) corresponds to a hydrodynamic limit. For the two different limits, several contributions have been performed. For instance, uniformly accurate numerical schemes have been proposed for the finite Larmor radius limit in [8], using two-scale strategies developed

<footnote>*Nantes Université, CNRS, Laboratoire de Mathématiques Jean Leray, LMJL, UMR 6629, F-44000 Nantes, France & Inria (MIN-GuS), France</footnote>

<footnote>†Univ Rennes, Inria (MINGuS) & IRMAR UMR 6625 & ENS Rennes, France</footnote>
in [5–7, 9, 14, 15]. On the other side, Asymptotic Preserving numerical schemes efficient in the hydrodynamic limit have been widely studied [10, 11, 17, 19, 24, 25, 27]. However, when two limits are combined, the literature is less abundant (see [13, 28]).

To investigate these stiffnesses, from a numerical point of view, a Particle-In-Cell (PIC) method is employed for the six-dimensional phase-space semi-discretization of (1). This method (see [1]) assumes the distribution function is approximated by a sum of macro-particles whose positions $x(t)$, velocities $v(t)$, and weights $\omega(t)$ satisfy a set of ODEs which, in our case, turns out to be a multiscale (highly-oscillatory and highly dissipative) set of ordinary differential equations (ODEs). Our goal is then to design efficient numerical time integrators to handle these multiscale characters using PIC method. On the one side, observing that particles positions and velocities are highly-oscillatory in time only, uniformly accurate (UA) time integrators can be derived thanks to the two-scale framework [5, 9, 14]. On the other side, the particles weights, which traduce the collision part, are solution to a nonlinear equation which combines the two stiff behaviors, highly oscillatory and dissipative. Roughly speaking, considering a BGK collision operator for $Q$ in (1), the particles weights $\omega(t)$ are solution to an equation of the form $\dot{\omega}(t) = \frac{1}{\tau}(M(t/\epsilon, \omega(t)) - \omega(t))$ with $M(\theta, \cdot)$ a nonlinear function of $\omega$. In this context, different multiscale numerical schemes are designed and analysed.

We need to consider the two-scale form of the weights since the highly-oscillatory behavior of the particles is transferred to the weights. Then, one of the main difficulties lies in the fact that the equilibria of the averaged collision operator are modified by the averaging procedure and it is crucial to take this asymptotic property into account in the numerical schemes. Indeed, in spite of its apparent simple relaxation form, the equilibria of the averaged collision operator are the so-called gyromaxwellian and have to be included in the numerical scheme to ensure its good asymptotic behavior. Following the penalty method combined with exponential integrators [17, 26], it is possible to design a scheme which enjoys the uniformly accuracy property with respect to $\epsilon$ (oscillations) and the asymptotic preserving property with respect to $\tau$ (relaxation). Moreover, in a reduced case where the filtered Maxwellian does not depend on the weights, it is possible to prove that the numerical scheme enjoys the uniform accuracy in both $\epsilon$ and $\tau$.

The rest of the paper is organized as follows. First, the different asymptotic models obtained from (1) are presented in Section 2: the averaging procedure leads to the effective collisional gyrokinetic model and the fluid asymptotic of this model leads to a gyrofluid model. Then, the multiscale numerical schemes are presented in Section 3 and analysed in Section 4. Lastly in Section 5, numerical results are presented to illustrate the properties of the numerical schemes in different configurations.

# 2 Presentation of the models

In this section, the different asymptotic models of (1) are recalled from [4]. We first recall the highly oscillatory limit $\epsilon \to 0$ ($\tau$ fixed): after filtering out the main oscillation, a collisional gyrokinetic model is obtained by a standard averaging procedure. The main point highlighted in [4] lies in the fact that the averaged collision operator $\langle Q \rangle$ has modified properties compared to the original collision operator $Q$. This is recalled in a second part. This paves the way of a new gyrofluid model whose properties are studied in Appendix A.

## 2.1 Gyrokinetic model

To formally derive the collisional gyrokinetic model from the limit $\epsilon \to 0$ of (1), we make explicit the $\epsilon$ and $\tau$ dependence in the distribution function unknown by denoting it $f^{\epsilon, \tau}(t, z)$ with $z = (x, v)$, so that (1) is recast as

$$ \partial_t f^{\epsilon, \tau} + \frac{1}{\epsilon} Az \cdot \nabla_z f^{\epsilon, \tau} + h(t, z) \cdot \nabla_z f^{\epsilon, \tau} = \frac{1}{\tau} Q(f^{\epsilon, \tau}), \quad (2) $$

with

$$ A = \begin{pmatrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \qquad h(t, z) = \begin{pmatrix} 0 \\ 0 \\ v_3 \\ E_1(t, x) \\ E_2(t, x) \\ E_3(t, x) \end{pmatrix}. \quad (3) $$

On the one side, the perpendicular (to the magnetic field) dynamics $\frac{1}{\epsilon}A$ generates oscillations caused by purely imaginary eigenvalues of $A$. On the other side, the parallel dynamics $h$ depends on the electric field $E = (E_1, E_2, E_3)$.
The first step is to perform a change of variable in order to filter out the main oscillation

$$Z = e^{-\frac{t}{\epsilon}A}z, \qquad (4)$$

and denote by $X = (X_1, X_2, X_3)$ the first 3 coordinates of $Z$ and by $V = (V_1, V_2, V_3)$ the last 3 coordinates of $Z$. This filtration is related to the flow $\mathfrak{z}$ associated to $Az \cdot \nabla_z$ by the relation $\mathfrak{z}(s, Z) = e^{sA}Z$. The flow is given by

$$ \begin{aligned} \mathfrak{z} &= (\mathfrak{x}, \mathfrak{y}), \\ \mathfrak{x}(s, Z) &= (\bar{X} + (I_2 - \mathcal{R}(-s))^{\perp} \bar{V}, X_3), \\ \mathfrak{y}(s, Z) &= (\mathcal{R}(-s) \bar{V}, V_3), \end{aligned} \qquad (5) $$

where $\bar{X} = (X_1, X_2)$, $\bar{V} = (V_1, V_2)$, $^{\perp}\bar{V} = (V_2, -V_1)$, $I_2$ the 2x2 identity matrix and $\mathcal{R}(\theta)$ stands for the rotation of angle $\theta$. Due to the spectral properties of $A$, one can prove that the flow $s \to \mathfrak{z}(s, Z)$ is $2\pi$-periodic, which will have important consequences from a numerical point of view. Let now define $F^{\epsilon, \tau}$ the filtered distribution function $F^{\epsilon, \tau}(t, Z) = f^{\epsilon, \tau}(t, z)$ (with $z = e^{t/\epsilon}A Z$ from (4)) which satisfies the filtered version of (2)

$$ \partial_t F^{\epsilon, \tau} + H(t, t/\epsilon, Z) \cdot \nabla_Z F^{\epsilon, \tau} = \frac{1}{\tau} Q_{\text{filt}}[F^{\epsilon, \tau}](t, t/\epsilon, Z), \qquad (6) $$

where

$$ Q_{\text{filt}}[F^{\epsilon, \tau}](t, s, Z) = Q[F^{\epsilon, \tau}(t, e^{-sA}.)](e^{sA}Z), \qquad H(t, s, Z) = e^{-sA}h(t, e^{sA}Z). \qquad (7) $$

As one can see, the definition of $Q_{\text{filt}}$ requires two changes of variables. First, the inverse change of variables enables to compute the moments in the original variables and then, the operator is evaluated in the new variables. As we will see in the next section, this will induce difficulties in the numerical approximation of $Q_{\text{filt}}$.

The filtered model (6) highlights three time scales. First, the slow time scale corresponding to $t \to H(t, s, Z)$ since the parallel dynamics evolves on a larger time scale than the perpendicular one; second the fast cyclotronic periodic time scale corresponding to $\theta \to H(t, \theta/\epsilon, Z)$; lastly, the presence of the collision term induces a collisional time scale $t/\tau$. The averaged model obtained when $\epsilon$ goes to 0 is derived by averaging with respect to the fast time variable $t/\epsilon$. Following [4], the collisional gyrokinetic model satisfied by $F^\tau$ (which is the limit as $\epsilon \to 0$ of $F^{\epsilon, \tau}$) writes

$$ \partial_t F^\tau + \langle H \rangle \cdot \nabla_Z F^\tau = \frac{1}{\tau} \langle Q \rangle [F^\tau], \qquad (8) $$

with

$$ \langle H \rangle (t, Z) = \frac{1}{2\pi} \int_{0}^{2\pi} H(t, s, Z) \, ds, \qquad \langle Q \rangle [F^{\tau}](t, Z) = \frac{1}{2\pi} \int_{0}^{2\pi} Q_{\text{filt}}[F^{\tau}](t, s, Z) \, ds. \qquad (9) $$

As (8) is still a collisional Vlasov equation, the hydrodynamic limit can be investigated and a fluid model can be derived. To do so, a careful study (performed in [4]) of the equilibrium and collision invariants of $\langle Q \rangle$ is necessary since, as $\tau \to 0$, $F^\tau$ converges to an equilibrium of $\langle Q \rangle$.

## 2.2 The average operator $\langle Q_{\text{BGK}} \rangle$

In view of numerical simulations, we will use the standard BGK operator which is the easiest operator satisfying the properties required by Theorem 4.1 of [4]. Indeed, the analysis presented in [4] does not depend on the collision operator choice, as long as it satisfies some properties. It is classical to check that the following BGK operator satisfies these properties (H-theorem, Maxwellian equilibria, moments conservation). In this subsection, we omit the time dependence for simplicity. First, we recall definition of the BGK operator:

$$ Q_{\text{BGK}}(f) = \mathcal{M}[f] - f, \qquad (10) $$

where the maxwellian $\mathcal{M}[f]$ associated to $f$ is defined by

$$ \mathcal{M}[f](x, v) = \frac{n(x)}{(2\pi T(x))^{3/2}} \exp\left(-\frac{|v - u(x)|^2}{2T(x)}\right), \qquad (11) $$
where the moments $n$, $u$ and $T$ are defined by

$$ 
\begin{aligned}
n(x) &= \int_{\mathbb{R}^3} f(x, v) \, dv, \qquad nu(x) = \int_{\mathbb{R}^3} f(x, v) v \, dv, \\
n\left(\frac{|u|^2}{2} + \frac{3}{2}T\right)(x) &= \int_{\mathbb{R}^3} f(x, v) \frac{|v|^2}{2} \, dv.
\end{aligned}
$$

In this section, we make the filtered model (6) and the averaged model (8) more explicit under the choice of $Q_{BGK}$. First, we define the filtered Maxwellian associated to a filtered density $F = F(Z)$ ($Z$ being defined in (4))

$$ \mathcal{M}_{\text{filt}}[F](s, Z) = \mathcal{M}[F(e^{-sA}.)](e^{sA}Z), \qquad (12) $$

and the averaged Maxwellian by

$$ \langle \mathcal{M} \rangle [F](Z) = \frac{1}{2\pi} \int_{0}^{2\pi} \mathcal{M}_{\text{filt}}[F](s, Z) \, ds. \qquad (13) $$

We give now the expression of the average operator $\langle Q_{BGK} \rangle$

$$ \langle Q_{\text{BGK}} \rangle [F](Z) = \langle \mathcal{M} \rangle [F](Z) - F(Z). \qquad (14) $$

Thanks to Proposition 4.5 of [4], the equilibria of this averaged BGK operator are the so-called gyromaxwellians defined by (recalling $\bar{X} = (X_1, X_2)$ and $\bar{V} = (V_1, V_2)$)

$$ \mathcal{G}[F](\bar{X}, X_3, V) = n \mathcal{M}_{\frac{\mu}{\mu-\theta}}^{2} (\bar{V} - \bar{U}) \mathcal{M}_{\theta}^{1}(V_3 - U_3) \mathcal{M}_{\mu}^{2}(\bar{X} + \perp \bar{V} - \bar{Y}), \qquad (15) $$

with $\mathcal{M}_T^d(v) = (2\pi T)^{-d/2} \exp(-|v|^2/2T)$ where $n(X_3) \in \mathbb{R}_+$, $U(X_3) \in \mathbb{R}^3$, $\bar{Y}(X_3) \in \mathbb{R}^2$, $\theta(X_3) \in \mathbb{R}^+$ and $\mu(X_3) \in \mathbb{R}_+$ are defined by integrating $F$ against the 8 collisional invariants of $\langle Q_{BGK} \rangle$ (whereas $Q_{BGK}$ has 5 collisional invariants)

$$ n(X_3) = \iint_{\mathbb{R}^2 \mathbb{R}^3} F(X, V) \, dV d\bar{X}, \qquad (16a) $$

$$ nU(X_3) = \iint_{\mathbb{R}^2 \mathbb{R}^3} F(X, V) V \, dV d\bar{X}, \quad n\bar{Y}(X_3) = \iint_{\mathbb{R}^2 \mathbb{R}^3} F(X, V)(\bar{X} + \perp \bar{V}) \, dV d\bar{X}, \qquad (16b) $$

$$ n\left(K + \frac{|U|^2}{2}\right)(X_3) = \iint_{\mathbb{R}^2 \mathbb{R}^3} F(X, V) \frac{|V|^2}{2} \, dV d\bar{X}, \qquad (16c) $$

$$ n\left(G + \frac{|\bar{Y}|^2}{2} - \frac{|\bar{U}|^2}{2}\right)(X_3) = \iint_{\mathbb{R}^2 \mathbb{R}^3} F(X, V) \frac{|\bar{X} + \perp \bar{V}|^2 - |\bar{V}|^2}{2} \, dV d\bar{X}. \qquad (16d) $$

These quantities, which only depend on $X_3$, are called gyromoments: $n$ is the density of particles, $U = (U_1, U_2, U_3)$ is the average velocity, $\bar{Y} = (Y_1, Y_2)$ the Larmor-center velocity, $K$ the internal kinetic energy and $G$ the internal gyration energy. Notice that $n > 0$, $K > 0$ and $K + G > 0$. The two temperatures $\theta$ and $\mu$ then verify $\mu > \theta > 0$ and satisfy the following relations

$$ K = \frac{\mu\theta}{\mu - \theta} + \frac{\theta}{2}, \qquad G = \mu - \frac{\mu\theta}{\mu - \theta}. \qquad (17) $$

We can inverse the system to express $\theta$ and $\mu$ as a function of $K$ and $G$ to get

$$ \theta = 2(K + G)/(1 + 2\nu), \qquad \mu = \nu\theta, \qquad (18) $$

with

$$ \nu = 1 + \frac{3S + \sqrt{9S^2 + 16(S+1)}}{4}, \qquad S = \frac{G}{K}. \qquad (19) $$

Moreover, we have the following relation: $\mu + \frac{\theta}{2} = K + G$.
**Remark 1.** We highlight the fact that $F \mapsto \langle Q_{BGK} \rangle [F] = \langle \mathcal{M} \rangle [F] - F$ and the relaxation operator $F \mapsto \mathcal{G}[F] - F$ are two different operators. As mentioned in [4], the average operator $\langle Q_{BGK} \rangle$ is more complicated than a relaxation type operator defined with gyromoments.

**Remark 2.** Let us remark that $\bar{P} := \bar{x} + \perp \bar{v}$ is invariant under the change of variable i.e. $\bar{x} + \perp \bar{v} = \bar{X} + \perp \bar{V}$. Thus, in this set of variables $(\bar{P}, V)$, we have $d\bar{x} \, d\bar{v} = d\bar{X} \, d\bar{V} = d\bar{P} \, d\bar{V}$ and the above calculations can be performed by replacing $F(X, V)$ by $\bar{F}(\bar{P}, X_3, V)$ in (16) and the gyromaxwellian can be written in a more natural way by

$$ \tilde{\mathcal{G}}[F](\bar{P}, X_3, V) = n \mathcal{M}_{\mu \theta}^{2} (\bar{V} - \bar{U}) \mathcal{M}_{\theta}^{1}(V_3 - U_3) \mathcal{M}_{\mu}^{2}(\bar{P} - \bar{Y}). \quad (20) $$

The variables $(\bar{P}, X_3, V)$ turn out to be natural for the gyrokinetic model (8).

A study of the gyrofluid model is presented in appendix A.

# 3 Discretisation of models

In this section, we develop multiscale numerical schemes for the kinetic model (1) or equivalently (2). The goal is to design numerical schemes which overcome the stiffness induced by the two physical parameters $\varepsilon$ and $\tau$. The main ingredients are the following: we discretize the phase space with a weighted particle method which transforms (2) into a set of ODEs; then, we construct uniformly accurate scheme in $\varepsilon$ using scales-separation strategy developed in [5–9, 14, 15] which introduces a new time periodic variable; the dissipative stiffness in $\tau$ is handled using IMEX and exponential time integrations which are widely used to achieve asymptotic preserving properties [17, 22, 26].

## 3.1 Particle-In-Cell framework

In order to discretize (2), we use a particle method. In the sequel, we will omit the notation referring to the solution dependence in $\tau$ and $\varepsilon$ to alleviate notations. Let $N_p \in \mathbb{N}$ a number of numerical macro-particles of positions $x_q$, velocities $v_q$ and weights $\omega_q$, $q = 1, \dots, N_p$. We will approximate the solution $f$ of (2) by the following particle ansatz

$$ f_{N_p}(t, x, v) = \sum_{q=1}^{N_p} \omega_q(t) \delta(x - x_q(t)) \delta(v - v_q(t)) = \sum_{q=1}^{N_p} \omega_q(t) \delta(z - z_q(t)), \quad (21) $$

with $z = (x, v)$ and $z_p(t) = (x_p(t), v_p(t))$. Evolution of particle $z_p(t)$ is given by the ODE system below, which corresponds to characteristics of (2).

$$ \dot{z}_p(t) = \frac{1}{\varepsilon}Az_p(t) + h(t, z_p(t)), \quad z_p(0) = z_{p,0}. \quad (22) $$

We assume that $z_{p,0} \in \Omega \cup D_v \subset \mathbb{R}^6$ where $\Omega \subset \mathbb{R}^3$ is the spatial domain and $D_v \subset \mathbb{R}^3$ the numerical velocity domain. Due to the presence of the source term $Q[f](t, z) = \frac{1}{\tau}(\mathcal{M}[f](t, z) - f(t, z))$, the weights $\omega_q(t)$ in (21) are not constant in time. We denote by $s_p(t) = \frac{1}{\tau}(m_p(t) - \omega_p(t))$ the weights of the source term with the Maxwellian's weights $m_p(t) = \mathcal{M}[f](t, z_p(t))V_{N_p}$ and $\omega_p(t) = f(t, z_p(t))V_{N_p}$ the weights of $f$, with $V_{N_p} = |\Omega||D_v|/N_p$ is the average volume of the phase space domain occupied by one particle. To compute the moments of $f_{N_p}$, we need to regularize $f_{N_p}$ with respect to $x$. In this work we consider tensorial product of first order B-splines: $S(x_1, x_2, x_3) = \tilde{S}(x_1)\tilde{S}(x_2)\tilde{S}(x_3)$. If $\Omega = [x_{min}, x_{max}]^3$ and $N_x$ the grid point number, we choose $\Delta x = (x_{max} - x_{min})/N_x$, then we have

$$ \tilde{S}(x) = \begin{cases} \frac{1}{\Delta x} \left( 1 - \frac{|x|}{\Delta x} \right) & \text{if } |x| \le \Delta x, \\ 0 & \text{else.} \end{cases} \quad (23) $$

With the smooth kernel $S$, we define the regularized form $f_{N_p}^S$ of $f_{N_p}$

$$ f_{N_p}^S(t, x, v) = \sum_{q=1}^{N_p} \omega_q(t) S(x - x_q(t)) \delta(v - v_q(t)), \quad (24) $$
from which we can compute the moments of $f_{N_p}^S$ following

$$ \left( n, nu, n \left( \frac{|u|^2}{2} + \frac{3}{2}T \right) \right) (t, x) = \sum_{q=1}^{N_p} \omega_q(t) S(x - x_q(t))(1, v_q(t), |v_q(t)|^2/2), \quad (25) $$

and then, we can define $\mathcal{M}[f_{N_p}^S]$ (using (11)) together with its weights

$$ m_p(t) = \mathcal{M}[f_{N_p}^S](t, z_p(t))V_{N_p}. \quad (26) $$

Notice that $m_p(t)$ does not only depend on $z_p$, but also depends on all weights and particles through $f_{N_p}^S$ and (24). To be consistent with the continuous case, we call the first one "functional" dependence (dependence with respect to $z_p$) and the second one "operational" dependence (dependence with respect to all $z_p$ and all $\omega_p$, which determine $f_{N_p}^S$). Finally, the weights $\omega_p$ satisfy the following equation

$$ \dot{\omega}_p(t) = \frac{1}{\tau}(m_p(t) - \omega_p(t)). \quad (27) $$

Now if we replace collisional operator by its particles version $\frac{1}{\tau}\sum_q(m_q - \omega_q)\delta(z - z_q)$ in (2), we can show, using (22)-(27), that particle density (21) is solution of (2) in the sense of distribution (see [12, 23, 31]). The goal is then to solve numerically the multiscale system (22)-(27) endowed with initial conditions $(z_p(0), \omega_p(0)) = (z_{p,0}, \omega_{p,0})$.

As in the continuous case, the main oscillation generated by $\dot{z}_p = (1/\epsilon)Az_p$ is filtered out by a change of unknown. We set $Z_p(t) = e^{-t/\epsilon}A z_p(t)$ (with $A$ given by (3)) so that (22) becomes

$$ \dot{Z}_p(t) = H(t, t/\epsilon, Z_p(t)) := e^{-t/\epsilon} h(t, e^{t/\epsilon} Z_p(t)). \quad (28) $$

Now, we explain how to compute $m_p(t)$ in (27) with filtered particules $Z_p(t) = (X_p(t), V_p(t))$. In this framework, the filtered distribution function writes

$$ F_{N_p}(t, Z = (X, V)) = \sum_{q=1}^{N_p} \omega_q(t) \delta(X - X_q(t)) \delta(V - V_q(t)). $$

To compute the moments, we need to apply the inverse change of variable. By using notations (5), we write the regularized form $F_{N_p}^S$ of $F_{N_p}$ thanks to the unfiltered particle density $f_{N_p}^S$ with filtered variables:

$$ \begin{aligned} F_{N_p}^S(t, Z) &= f_{N_p}^S(t, \mathfrak{X}(t/\epsilon, Z), \mathfrak{V}(t/\epsilon, Z)) \\ &= \sum_{q=1}^{N_p} \omega_q(t) S(\mathfrak{X}(t/\epsilon, Z) - \mathfrak{X}(t/\epsilon, Z_q(t))) \delta(\mathfrak{V}(t/\epsilon, Z) - \mathfrak{V}(t/\epsilon, Z_q(t))). \end{aligned} \quad (29) $$

We are now able to compute moments of $F_{N_p}^S$ as

$$ \begin{aligned} \mathcal{N}(t, Z) &= \sum_{q=1}^{N_p} \omega_q(t) S(\mathfrak{X}(t/\epsilon, Z) - \mathfrak{X}(t/\epsilon, Z_q(t))), \\ \mathcal{NU}(t, Z) &= \sum_{q=1}^{N_p} \omega_q(t) S(\mathfrak{X}(t/\epsilon, Z) - \mathfrak{X}(t/\epsilon, Z_q(t))) \mathfrak{V}(t/\epsilon, Z_q(t)), \end{aligned} \quad (30) $$

$$ \mathcal{N}\left(\frac{|\mathcal{U}|^2}{2} + \frac{3}{2}\mathcal{T}\right)(t, Z) = \sum_{q=1}^{N_p} \omega_q(t) S(\mathfrak{X}(t/\epsilon, Z) - \mathfrak{X}(t/\epsilon, Z_q(t))) \frac{|\mathfrak{V}(t/\epsilon, Z_q(t))|^2}{2}, $$

and to reconstruct the corresponding filtered Maxwellian $\mathcal{M}_{\text{filt}}[F_{N_p}^S]$. Let us remark that we have $(\mathcal{N}, \mathcal{U}, \mathcal{T})(t, Z) = (n, u, T)(t, \mathfrak{X}(t/\epsilon, Z))$ and that the weights of the Maxwellian $\mathcal{M}$ and the filtered Maxwellian $\mathcal{M}_{\text{filt}}[F]$ are equal since the filtration only applies on the particles

$$ \mathcal{M}[f_{N_p}^S](t, z) = \sum_{q=1}^{N_p} m_p(t) \delta(z - z_q(t)), \quad \mathcal{M}_{\text{filt}}[F_{N_p}^S](t, Z) = \sum_{q=1}^{N_p} m_p(t) \delta(Z - Z_q(t)), $$
and the same holds true for $f_{N_p}$ and $F_{N_p}$. Finally, once the filtered Maxwellian $\mathcal{M}_{\text{filt}}[F_{N_p}^S](t, Z)$ has been constructed, the weights are computed by evaluating it at the particles $Z_p(t)$: $m_p(t) = \mathcal{M}_{\text{filt}}[F_{N_p}^S](t, Z_p(t))V_{N_p}$. To anticipate with the next subsection devoted to the numerical schemes, we introduce a new definition to highlight the "operational" dependence of $\mathcal{M}_{\text{filt}}$ (with respect to the weights and particles positions) and its dependence with respect to the fast variable

$$m_p(t) = \mathcal{M}_{\text{filt}}[Z(t), \omega(t)](t, t/\epsilon, Z_p(t))V_{N_p}, \quad (31)$$

where $Z(t) = \{Z_q(t)\}_{q=1,...,N_p}$ and $\omega(t) = \{\omega_q(t)\}_{q=1,...,N_p}$.

## 3.2 Multiscale numerical schemes

The presence of a slow time scale and a fast one which is $2\pi$-periodic is amenable to the two-scale framework [5-9,14,15]. Then, we consider the two time scales as independent by introducing $\mathcal{Z}_p(t, s)$ (resp. $\mathcal{W}_p(t, s)$) satisfying $\mathcal{Z}_p(t, t/\epsilon) = Z_p(t)$ (resp. $\mathcal{W}_p(t, t/\epsilon) = \omega_p(t)$). We also introduce the following functions:

$$\begin{aligned} \mathcal{H}_p(t, s) &= H(t, s, \mathcal{Z}_p(t, s)), \\ \mathcal{M}_p(t, s) &= \mathcal{M}_{\text{filt}}[\mathcal{Z}(t, s), \mathcal{W}(t, s)](t, s, \mathcal{Z}_p(t, s))V_{N_p}. \end{aligned} \quad (32)$$

With these notations, (28) and (27) become

$$\begin{aligned} \partial_t \mathcal{Z}_p(t, s) + \frac{1}{\epsilon} \partial_s \mathcal{Z}_p(t, s) &= \mathcal{H}_p(t, s), \\ \partial_t \mathcal{W}_p(t, s) + \frac{1}{\epsilon} \partial_s \mathcal{W}_p(t, s) &= \frac{1}{\tau} (\mathcal{M}_p(t, s) - \mathcal{W}_p(t, s)), \end{aligned} \quad (33)$$

where the calculation of the moments, the Maxwellian and its weights presented in Subsection 3.1 are adapted to the two-scale framework.

Let us now discuss the time integration of (33). In this part, we present different multiscale numerical schemes according to the considered asymptotic in $\epsilon$ or/and $\tau$. Since quantities are periodic in $s$, we will consider the Fourier transform to approximate the functions in this direction. With a uniform mesh $s_j = j2\pi/N_s, j = 0, ..., N_s$, with $N_s \in \mathbb{N}^*$ (an even number), the Fourier transform writes, for an arbitrary function $g: s \in [0, 2\pi] \to g(s)$ of Fourier coefficients $\hat{g}_j$,

$$g(s) = \sum_{l=-N_s/2}^{N_s/2-1} \hat{g}_l e^{ils}.$$

Applying the Fourier transform to (33) enables to get the following equations on the modes

$$\begin{aligned} \frac{d}{dt} \hat{Z}_{p,l}(t) + \frac{il}{\epsilon} \hat{Z}_{p,l}(t) &= \hat{\mathcal{H}}_{p,l}(t), \\ \frac{d}{dt} \hat{W}_{p,l}(t) + \left( \frac{il}{\epsilon} + \frac{1}{\tau} \right) \hat{W}_{p,l}(t) &= \frac{1}{\tau} \hat{\mathcal{M}}_{p,l}(t). \end{aligned}$$

We introduce time discretisation $t^n = n\Delta t, \forall n \in \mathbb{N}$ with $\Delta t > 0$ the time step. Applying the variation of constant formula between $t^n$ and $t^{n+1}$ leads to

$$\begin{aligned} \hat{Z}_{p,l}(t^{n+1}) &= e^{-\frac{il}{\epsilon}\Delta t} \hat{Z}_{p,l}(t^n) + \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\epsilon}(t^{n+1}-t)} \hat{\mathcal{H}}_{p,l}(t) \mathrm{d}t, \\ \hat{W}_{p,l}(t^{n+1}) &= e^{-\left(\frac{il}{\epsilon} + \frac{1}{\tau}\right)\Delta t} \hat{W}_{p,l}(t^n) + \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-\left(\frac{il}{\epsilon} + \frac{1}{\tau}\right)(t^{n+1}-t)} \hat{\mathcal{M}}_{p,l}(t) \mathrm{d}t. \end{aligned} \quad (34)$$

These expressions will pave the way of the time integrators that will be described in the sequel.
### 3.3 Scheme A: uniform accuracy with respect to $\epsilon$

In this subsection, we introduce a first numerical scheme to approximate (33). The use of the two-scale approach will make this scheme uniformly accurate with respect to $\epsilon$ for all fixed $\tau > 0$. Besides, this scheme enjoys the asymptotic preserving property in the asymptotic $\tau \to 0$ for all fixed $\epsilon > 0$.

The starting point for the numerical scheme is (34). To discretize the integrals in (34), we use an Adam-Bashforth type method. We consider here first and second order methods but higher order can be considered as in [9], [29]. Denoting by $g(t)$ a function representing $\hat{\mathcal{H}}_{p,l}(t)$ or $\hat{\mathcal{M}}_{p,l}(t)$ in (34), we replace $g(t)$ by the following approximations

$$g(t) \approx g(t^n) \quad (\text{first order approximation}),$$

$$g(t) \approx g(t^n) + (t - t^n) \frac{g(t^n) - g(t^{n-1})}{\Delta t} \quad (\text{second order approximation}).$$

To simplify the final expression of the scheme, we introduce for all $\alpha \in \mathbb{N}$, $z \in \mathbb{C}$, the following function

$$\varphi_{\alpha}(z) := \frac{1}{\alpha!} \int_{0}^{1} e^{-z(1-t)} t^{\alpha} \, dt = \frac{(-1)^{\alpha+1}}{z^{\alpha+1}} \left( e^{-z} - \sum_{k=0}^{\alpha} \frac{(-z)^k}{k!} \right),$$

so that

$$\int_{t^n}^{t^{n+1}} e^{-z(t^{n+1}-t)} (t-t^n)^{\alpha} dt = \alpha! \Delta t^{\alpha+1} \varphi_{\alpha}(z\Delta t). \quad (35)$$

We denote $\mathcal{Z}_p^n(s)$ (resp. $\mathcal{W}_p^n(s)$) an approximation of $\mathcal{Z}_p(t^n, s)$ (resp. $\mathcal{W}_p(t^n, s)$), and $\hat{\mathcal{Z}}_{p,l}^n$ (resp. $\hat{\mathcal{W}}_{p,l}^n$) the coefficients of the Fourier transform of $\mathcal{Z}_p^n(s)$ (resp. $\mathcal{W}_p^n(s)$). With these above notations, the first order numerical scheme writes

$$\hat{\mathcal{Z}}_{p,l}^{n+1} = e^{-\frac{il}{\epsilon}\Delta t} \hat{\mathcal{Z}}_{p,l}^{n} + \Delta t \varphi_0 \left(\frac{il}{\epsilon}\Delta t\right) \hat{\mathcal{H}}_{p,l}^{n}, \quad (36a)$$

$$\hat{\mathcal{W}}_{p,l}^{n+1} = e^{-(\frac{il}{\epsilon}+\frac{1}{\tau})\Delta t} \hat{\mathcal{W}}_{p,l}^{n} + \frac{\Delta t}{\tau} \varphi_0 \left( \left( \frac{il}{\epsilon} + \frac{1}{\tau} \right) \Delta t \right) \hat{\mathcal{M}}_{p,l}^{n}, \quad (36b)$$

and the second order one is

$$\begin{aligned} \hat{\mathcal{Z}}_{p,l}^{n+1} &= e^{-\frac{il}{\epsilon}\Delta t} \hat{\mathcal{Z}}_{p,l}^{n} + \Delta t \varphi_0 \left( \Delta t \frac{il}{\epsilon} \right) \hat{\mathcal{H}}_{p,l}^{n} + \Delta t \varphi_1 \left( \Delta t \frac{il}{\epsilon} \right) (\hat{\mathcal{H}}_{p,l}^{n} - \hat{\mathcal{H}}_{p,l}^{n-1}), \\ \hat{\mathcal{W}}_{p,l}^{n+1} &= e^{-(\frac{il}{\epsilon}+\frac{1}{\tau})\Delta t} \hat{\mathcal{W}}_{p,l}^{n} + \frac{\Delta t}{\tau} \varphi_0 \left( \Delta t \left( \frac{il}{\epsilon} + \frac{1}{\tau} \right) \right) \hat{\mathcal{M}}_{p,l}^{n} \\ &\quad + \frac{\Delta t}{\tau} \varphi_1 \left( \Delta t \left( \frac{il}{\epsilon} + \frac{1}{\tau} \right) \right) (\hat{\mathcal{M}}_{p,l}^{n} - \hat{\mathcal{M}}_{p,l}^{n-1}). \end{aligned} \quad (37)$$

Then, once the two-scale (Fourier) numerical solutions are computed, one can evaluate the Fourier series

$$\mathcal{Z}_p^n(s) = \sum_{l=-N_s/2}^{N_s/2-1} \hat{\mathcal{Z}}_{p,l}^n e^{ils}, \qquad \mathcal{W}_p^n(s) = \sum_{l=-N_s/2}^{N_s/2-1} \hat{\mathcal{W}}_{p,l}^n e^{ils}, \quad (38)$$

and the original solutions are obtained by evaluating on the diagonal $s = t^n/\epsilon$ and applying the inverse change of variable

$$Z_p^n = Z_p^n \left( \frac{t^n}{\epsilon} \right), \qquad (x_p^n, v_p^n) = z_p^n = e^{\frac{t^n}{\epsilon} A} Z_p^n, \qquad \omega_p^n = W_p^n \left( \frac{t^n}{\epsilon} \right). \quad (39)$$

Let us discuss the initial condition. By introducing the variable $s$, we obtain a degree of freedom on the initial condition for (33). The only constraint is $\mathcal{Z}_p(0, 0) = Z_p(0)$ and $\mathcal{W}_p(0, 0) = \omega_p(0)$. We use this degree of freedom to control the smoothness of $\mathcal{Z}_p(t)$ and $\mathcal{W}_p(t)$. To do so, well-prepared initial conditions can be constructed that allow the scheme to be uniformly accurate with respect to $\epsilon$. For the first order approximation, we consider constant in $s$ initial condition

$$\mathcal{Z}_p(0, s) = Z_p(0) = (x_p(0), v_p(0)), \qquad \mathcal{W}_p(0, s) = \omega_p(0). \quad (40)$$
For the second order approximation, we consider

$$
\begin{align}
\mathcal{Z}_p(0, s) &= Z_{p,0} + \varepsilon \int_0^s (I - \Pi) H(0, \theta, Z_{p,0}) \, d\theta, \tag{41} \\
\mathcal{W}_p(0, s) &= \omega_{p,0} + \frac{\varepsilon}{\tau} \int_0^s (I - \Pi) \mathcal{M}_p(0, \theta) \, d\theta,
\end{align}
$$

where $\Pi g = \frac{1}{2\pi} \int_0^{2\pi} g(s) ds$ is the averaging operator (see [5] for more details). This preparation is only required to ensure uniform estimates of the time derivatives with respect to $\varepsilon$ and is not used in the classical fluid limit ($\varepsilon > 0$ fixed and $\tau \to 0$). The two-scale framework provides uniform accuracy for (36) and (37) with respect to $\varepsilon$, for fixed $\tau$. Conversely, for a fixed $\varepsilon$ and looking for the limit $\tau \to 0$, the problem becomes nonstiff in $\varepsilon$ and dissipative in $\tau$. For such a problem, several numerical schemes have been proposed, entering in the category of Asymptotic Preserving schemes [24]. This is the purpose of the next part.

3.4 Scheme A: Asymptotic property with respect to τ

In this subsection, the asymptotic behavior of the first order scheme (36) is investigated. Let us recall the scheme
on the weights (36b) (the equation on the particles (36a) being unchanged)

$$
\hat{W}_{p,l}^{n+1} = e^{-(\frac{il}{\varepsilon} + \frac{1}{\tau})\Delta t} \hat{W}_{p,l}^{n} + \frac{\varepsilon}{il\tau + \varepsilon} \left(1 - e^{-(\frac{il}{\varepsilon} + \frac{1}{\tau})\Delta t}\right) \hat{M}_{p,l}^{*}
$$

with $\hat{\mathcal{M}}_{p,l}^* = \hat{\mathcal{M}}_{p,l}^n$. We have an AP property in the following sense: when $\tau \to 0$, we have $\hat{W}_{p,l}^{n+1} \to \hat{\mathcal{M}}_{p,l}^*$ and thus our density weights converge toward maxwellian weights. However to get a better asymptotic behavior when $\tau \to 0$, an implicit scheme is required by taking $\hat{\mathcal{M}}_{p,l}^* = \hat{\mathcal{M}}_{p,l}^{n+1}$. However, computing $\hat{\mathcal{M}}_{p,l}^{n+1}$ requires the knowledge of $\hat{Z}_{p,l}^{n+1}$ and $\hat{W}_{p,l}^{n+1}$ according to (31). Actually, $\hat{\mathcal{M}}_{p,l}$ depends on the weights through the moments, but here, contrary to the usual Maxwellian case [11,24], the moments (30) are not constant during this (collisional) step. Indeed, to compute the filtered Maxwellian, the moments of (2) have to be considered and we observe that not only they evolve in time but they are stiff with respect to $\varepsilon$. Thus, to avoid a costly fully implicit scheme, we consider a hybrid approximation of (31) computed from the knowledge of $\hat{Z}_{p,l}^{n+1}$ and $\hat{W}_{p,l}^n$.

The procedure is then the following: we first compute $\mathcal{Z}_q^{n+1}(s)$ from the first equation of (36) and then we
compute the moments according to

$$
\mathbf{U}^*(s, x) = \sum_{q=1}^{N_p} \mathcal{W}_p^n(s) S\left(x - \mathcal{X}\left(s + \frac{\Delta t}{\varepsilon}, \mathcal{Z}_q^{n+1}\left(s + \frac{\Delta t}{\varepsilon}\right)\right)\right) \bar{m}\left(\mathfrak{V}\left(s + \frac{\Delta t}{\varepsilon}, \mathcal{Z}_q^{n+1}\left(s + \frac{\Delta t}{\varepsilon}\right)\right)\right), \quad (42)
$$

where $\bar{m}(v) = (1, v, |v|^2/2)$. Then, we are able to compute the Maxwellian weights $\mathcal{M}_p^*(s)$ according to these
moments:

$$
\mathcal{M}_{p}^{*}(s) = \mathcal{M}_{\text{filt}}[\mathcal{Z}^{n+1}(s + \Delta t/\epsilon), \mathcal{W}^{n}(s)](t^{n+1}, s + \Delta t/\epsilon, \mathcal{Z}_{p}^{n+1}(s + \Delta t/\epsilon))V_{N_{p}}. \quad (43)
$$

Once the weights of the filtered Maxwellian are computed, we then update the weights with

$$
\hat{W}_{p,l}^{n+1} = e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t} \hat{W}_{p,l}^{n} + \frac{\epsilon}{il\tau + \epsilon} \left(1 - e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t}\right) \hat{M}_{p,l}^{*} \quad (44)
$$

with $\hat{\mathcal{M}}_{p,l}^*$ computed with (43). Then, when $\tau \to 0$, $\hat{W}_{p,l}^{n+1}$ converges towards $\hat{\mathcal{M}}_{p,l}^*$.

For the second order scheme, when $\tau$ goes to 0, $\hat{W}_{p,l}^{n+1}$ converges toward $2\hat{\mathcal{M}}_{p,l}^n - \hat{\mathcal{M}}_{p,l}^{n-1}$, which is the linear extrapolation of $\hat{\mathcal{M}}_{p,l}^{n+1}$.

3.5 Scheme B: AP properties in the gyrofluid limit

In the gyrofluid limit, the distribution function converges towards the gyromaxwellian (15) which belongs to the
kernel of the averaged BGK operator $\langle Q_{\text{BGK}} \rangle$. This important property is not ensured by the scheme A presented
in Subsections 3.3 and 3.4. When $\epsilon \to 0$, the averaged BGK operator $\langle Q_{BGK} \rangle [F] = \langle M_{\text{filt}} \rangle [F] - F$ turns out to be more complicated than standard BGK-like operator. Indeed, as $\tau \to 0$, $F$ converges towards $\mathcal{G}[F]$ and not $\langle M_{\text{filt}} \rangle [F]$. Hence, the averaged BGK operator $\langle Q_{BGK} \rangle$ shares similarities with more complex nonlinear operators like the Boltzmann or Landau operator for which specific Asymptotic Preserving schemes have been developed (see [17, 18]) to ensure the relaxation of the distribution function towards the Maxwellian equilibria which are not captured by standard schemes. A similar strategy is used here to enforce the relaxation of the scheme to the gyromaxwellian $\mathcal{G}[F]$ defined in (15).

In this section, we focus on the relaxation part ($\tau \to 0$) and thus, we only consider the equation on the weights, the particles positions being approximated by (36a). Starting from the equation (33) on the weights, we add and subtract the gyromaxwellian weights (denoted by $(1/\tau) \mathcal{G}_p$). We compute $\mathcal{G}_p(t, s)$ in a similar way than maxwellian weights: we compute gyromoments associated to $\mathcal{W}_p$, which determines a unique gyromaxwellian $\mathcal{G}$. Then we compute gyromaxwellian weights as in (32) with $\mathcal{G}$ instead of $M_{\text{filt}}$. Equation (33) on the weights becomes

$$ \partial_t \mathcal{W}_p + \frac{1}{\epsilon} \partial_s \mathcal{W}_p = \frac{1}{\tau}(M_p - \mathcal{W}_p) - \frac{1}{\tau} \mathcal{G}_p + \frac{1}{\tau} \mathcal{G}_p. \quad (45) $$

Considering the Fourier transform in $s$, the variation of constant formula leads to

$$ \partial_t (\hat{\mathcal{W}}_{p,l} e^{(\frac{il}{\epsilon} + \frac{1}{\tau})t}) = e^{(\frac{il}{\epsilon} + \frac{1}{\tau})t} \left[ \frac{1}{\tau}(\hat{M}_{p,l} - \hat{\mathcal{G}}_{p,l}) + \frac{1}{\tau}\hat{\mathcal{G}}_{p,l} \right]. \quad (46) $$

We now integrate between $t^n$ and $t^{n+1}$ to get

$$ \begin{aligned} \hat{\mathcal{W}}_{p,l}(t^{n+1}) = & e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t} \hat{\mathcal{W}}_{p,l}(t^n) + \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})(t^{n+1}-t)} (\hat{M}_{p,l}(t) - \hat{\mathcal{G}}_{p,l}(t)) dt \\ & + \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})(t^{n+1}-t)} \hat{\mathcal{G}}_{p,l}(t) dt, \end{aligned} $$

and perform different quadratures (Lawson and exponential integrators type strategies) on the two last right hand side terms as follows

$$ \begin{aligned} \hat{\mathcal{W}}_{p,l}^{n+1} = & e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t} \hat{\mathcal{W}}_{p,l}^n + \frac{\epsilon}{il\tau} e^{-\frac{\Delta t}{\tau}} (1 - e^{-\frac{il}{\epsilon}\Delta t}) (\hat{M}_{p,l}^n - \hat{\mathcal{G}}_{p,l}^n) \\ & + \frac{\epsilon}{il\tau + \epsilon} (1 - e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t}) \hat{\mathcal{G}}_{p,l}^{n+1}, \quad l \neq 0, \end{aligned} \quad (47) $$

$$ \hat{\mathcal{W}}_{p,0}^{n+1} = e^{-\frac{\Delta t}{\tau}} \hat{\mathcal{W}}_{p,0}^n + \frac{\Delta t}{\tau} e^{-\frac{\Delta t}{\tau}} (\hat{M}_{p,0}^n - \hat{\mathcal{G}}_{p,0}^n) + (1 - e^{-\frac{\Delta t}{\tau}}) \hat{\mathcal{G}}_{p,0}^{n+1}. $$

To compute $\hat{\mathcal{G}}_{p,l}^{n+1}$, we need to compute the gyromoments $\mathbf{U}^{n+1}$ defined by (16). To do so, we integrate (45) with respect to the invariants (assuming that particles have already evolved and are frozen during this step) to get

$$ \partial_t \mathbf{U} + \frac{1}{\epsilon} \partial_s \mathbf{U} = \frac{1}{\tau}(\tilde{\mathbf{U}} - \mathbf{U}), $$

where $\tilde{\mathbf{U}}$ denotes the gyromoments of the filtered Maxwellian $M_{\text{filt}}$ (whose weights are denoted by $M_p$). An exponential scheme to update the Fourier modes of $\mathbf{U}$ from $t^n$ to $t^{n+1}$ is given by

$$ \begin{aligned} \hat{\mathbf{U}}_{l}^{n+1} = & e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t} \hat{\mathbf{U}}_{l}^{n} + \frac{\epsilon}{il\tau + \epsilon} (1 - e^{-(\frac{il}{\epsilon} + \frac{1}{\tau})\Delta t}) \hat{\mathbf{U}}_{l}^{n}, \quad l \neq 0, \\ \hat{\mathbf{U}}_{0}^{n+1} = & e^{-\frac{\Delta t}{\tau}} \hat{\mathbf{U}}_{0}^{n} + (1 - e^{-\frac{\Delta t}{\tau}}) \hat{\mathbf{U}}_{0}^{n}. \end{aligned} \quad (48) $$

As we will see in the next section, this scheme is AP in the gyrokinetic limit ($\epsilon \to 0$ and $\tau > 0$ fixed) and in the gyrofluid limit ($\epsilon \to 0$ then $\tau \to 0$). However, for this scheme, the classical fluid limit ($\tau \to 0$ and $\epsilon > 0$ fixed) is not recovered since $\forall \epsilon \ge 0$, $\hat{\mathcal{W}}_{p,l}^{n+1} \to \hat{\mathcal{G}}_{p,l}^{n+1}$ when $\tau \to 0$.

## 3.6 Convex combined scheme

The scheme B presented in the last subsection does not capture the correct limit when $\tau \to 0$ for a fixed $\epsilon$ since it converges to a gyromaxwellian by construction. With the notation $\hat{\mathcal{W}}_{p,l}^{n+1,A}$ obtained from (44) and $\hat{\mathcal{W}}_{p,l}^{n+1,B}$
obtained from (47)-(48), we consider a convex combination of these two numerical schemes with coefficients $\frac{\epsilon}{\tau+\epsilon}$ and $(1 - \frac{\epsilon}{\tau+\epsilon})$ so that, according to the desired regimes, the coefficients degenerate towards the correct scheme:

$$ \hat{W}_{p,l}^{n+1} = \frac{\epsilon}{\tau + \epsilon} \hat{W}_{p,l}^{n+1,A} + \left(1 - \frac{\epsilon}{\tau + \epsilon}\right) \hat{W}_{p,l}^{n+1,B}. $$

To summarise the different asymptotics of this convex combined scheme, when $\epsilon > 0$ fixed and $\tau \to 0$, the scheme degenerates into scheme A which has right behaviour in the classical fluid limit. When $\tau > \epsilon \to 0$, the scheme degenerates into scheme B which has the good behaviour in the gyrofluid limit.

# 4 Asymptotic properties of the schemes

In this section, we study the asymptotic properties of the schemes presented above. As we will see in the first subsection, scheme B is AP in the gyrofluid limit ($\epsilon \to 0$ then $\tau \to 0$). In the second subsection, we show that, when the Maxwellian is independent of weights, scheme A is UA in $\epsilon$ and in $\tau$.

## 4.1 AP properties of scheme B for the gyrofluid limit

We investigate the gyrofluid limit of scheme B. The next proposition shows that scheme B enjoys the Asymptotic Preserving property in the gyrofluid limit.

**Proposition 1.** The scheme B defined by (36a)-(47)-(48) enjoys the Asymptotic Preserving property in the gyrofluid limit $\epsilon \to 0$ then $\tau \to 0$.

*Proof.* We first consider the limit $\epsilon \to 0$ in the scheme B, which corresponds to looking at the mode $l = 0$ in Fourier (observing that all the Fourier modes converge towards 0 as $\epsilon \to 0$ except for the mode $l = 0$). Then, in (36a), we obtain when $\epsilon \to 0$

$$ \hat{Z}_{p,0}^{n+1} = \hat{Z}_{p,0}^{n} + \Delta t \hat{H}_{p,0}^{n}, $$

whereas for the weights equations (47)-(48), we obtain when $\epsilon \to 0$

$$ \hat{W}_{p,0}^{n+1} = e^{-\frac{\Delta t}{\tau}} \hat{W}_{p,0}^{n} + \frac{\Delta t}{\tau} e^{-\frac{\Delta t}{\tau}} (\hat{M}_{p,0}^{n} - \hat{G}_{p,0}^{n}) + (1 - e^{-\frac{\Delta t}{\tau}}) \hat{G}_{p,0}^{n+1}, \quad (49) $$

$$ \hat{U}_{0}^{n+1} = e^{-\frac{\Delta t}{\tau}} \hat{U}_{0}^{n} + (1 - e^{-\frac{\Delta t}{\tau}}) \hat{U}_{0}^{n} = \hat{U}_{0}^{n}, \quad (50) $$

since $\hat{U}_{0}^{n} = \hat{U}_{0}^{n}$ when $\epsilon = 0$ (indeed, $\langle Q_{BGK} \rangle$ conserves gyromoments). This scheme is a consistent explicit first order time discretization of the gyrokinetic equation (8).

We now focus on the limit $\tau \to 0$. Following [11, 30, 33] in which an asymptotic preserving particle method is used to deal with the hydrodynamic limit, we check that the weights of the unknown degenerate when $\tau \to 0$ towards the gyromaxwellian weights. It is clear from (49)-(50) that $\hat{W}_{p,0}^{n+1}$ converges to $\hat{G}_{p,0}^{n+1}$ which means that the asymptotic scheme is a consistent discretization of the gyrofluid equations (65) (presented in Appendix A). $\square$

## 4.2 Uniform Accuracy for scheme A

In this section, we study the convergence of scheme A presented above. As the full collisional context remains difficult to analyse, we first consider the case where the Maxwellian is a given function depending on $\mathcal{Z}$ only. In this case, one can prove uniform accuracy in time with respect to both $\epsilon$ and $\tau$. When the Maxwellian depends on the weights $\mathcal{W}$, we are able to prove uniform accuracy in time with respect to $\epsilon$ (using classical tools from [5], [32]) and an asymptotic property in $\tau$.

The equation on $\mathcal{Z}_p$ enters in the two-scale framework proposed and analysed in [5–9, 14, 15]. The additional variable $s \in [0, 2\pi]$ offers a degree of freedom to choose the initial condition $\mathcal{Z}_p(0, s)$ which is well-prepared, i.e. $\mathcal{Z}_p(t, s)$ becomes smooth with respect to $\epsilon$. To be more precise, following a Chapman-Enskog type expansion, the first order initial condition chosen as (40) ensures that the time derivatives of $\mathcal{Z}_p(t, s)$ are uniformly bounded

$$ \sup_{t \in [0, t_{\text{max}}]} \| \partial_t^k \mathcal{Z}_p(t, \cdot) \|_{L^\infty} \le C, \quad k = 0, 1, 2. \quad (51) $$
We consider the model for $(\mathcal{Z}_p, \mathcal{W}_p) : (t, s) \in [0, t_{\max}] \times [0, 2\pi] \rightarrow (\mathcal{Z}, \mathcal{W})(t, s) \in \mathbb{R}^{6N_p} \times \mathbb{R}^{N_p}$

$$ 
\begin{aligned}
\partial_t \mathcal{Z}_p(t, s) + \frac{1}{\varepsilon} \partial_s \mathcal{Z}_p(t, s) &= \mathcal{H}_p(t, s), \\
\partial_t \mathcal{W}_p(t, s) + \frac{1}{\varepsilon} \partial_s \mathcal{W}_p(t, s) &= \frac{1}{\tau} (\mathcal{M}_p(t, s) - \mathcal{W}_p(t, s)),
\end{aligned}
\qquad (52) $$

where $\mathcal{H}_p(t, s) = H(t, s, \mathcal{Z}_p(t, s)) \in \mathbb{R}^{6N_p}$ and $\mathcal{M}_p(t, s) = M(t, s, \mathcal{Z}_p(t, s)) \in \mathbb{R}^{N_p}$, with a well-prepared initial condition $\mathcal{Z}_p(0, s), \mathcal{W}_p(0, s)$ computed as in (41). Obviously, some assumptions are required. In the sequel, the norm $\| \cdot \|$ of a vector refers to the standard euclidean norm in $\mathbb{R}^d$ ($d \in \mathbb{N}^*$). Some smoothness is required for $H$ and $M$ ($C^2$ regularity). Under these assumptions, it can be proved that (51) is true. Let us remark that we assumed here that $M$ is independent of $\omega$. This point will be discussed below.

**Theorem 1.** Under the above assumptions, the scheme (37)-(38)-(39)

with initial condition (41) is uniformly accurate in $\varepsilon$ and $\tau$: there exists a constant $C > 0$, independent of $\varepsilon$ and $\tau$ such that

$$ \|\mathcal{Z}^n(t^n/\varepsilon) - z(t^n)\|_{L^\infty} + \|\mathcal{W}^n(t^n/\varepsilon) - w(t^n)\|_{L^\infty} \le C\Delta t^2, \quad 0 \le n \le N, \quad N\Delta t = t_{\max}, \qquad (53) $$

where $z(t), w(t)$ are the solutions of (22) and (27).

*Proof.* Recall that

$$ 
\begin{aligned}
\mathcal{H}_p(t, s) &= H(t, s, \mathcal{Z}_p(t, s)), \\
\mathcal{M}_p(t, s) &= M_{\text{filt}}[\mathcal{Z}(t, s), \mathcal{W}(t, s)](t, s, \mathcal{Z}_p(t, s)),
\end{aligned}
 $$

and

$$ 
\begin{aligned}
\mathcal{H}_p^n(s) &= H(t^n, s, \mathcal{Z}_p^n(s)), \\
\mathcal{M}_p^n(s) &= M_{\text{filt}}[\mathcal{Z}^n(s), \mathcal{W}^n(s)](t^n, s, \mathcal{Z}_p^n(s)).
\end{aligned}
 $$

The solution $\mathcal{Z}_p$ of the two-scale equation (52) satisfies

$$ 
\begin{aligned}
\mathcal{Z}_p(t^{n+1}, s) &= \sum_{l \in \mathbb{Z}} \hat{\mathcal{Z}}_{p,l}(t^{n+1}) e^{ils} \\
&= \sum_{l \in \mathbb{Z}} \left( e^{-\frac{il}{\varepsilon} \Delta t} \hat{\mathcal{Z}}_{p,l}(t^n) + \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\varepsilon} (t^{n+1}-t)} \hat{\mathcal{H}}_{p,l}(t) \, dt \right) e^{ils}. \qquad (54)
\end{aligned}
 $$

A second order Taylor-Lagrange formula enables to expand the term $\hat{\mathcal{H}}_{p,l}(t)$ at $t = t^n$ to get

$$ 
\begin{aligned}
\mathcal{Z}_p(t^{n+1}, s) = & \sum_{l \in \mathbb{Z}} \left( e^{-\frac{il}{\varepsilon} \Delta t} \hat{\mathcal{Z}}_{p,l}(t^n) + \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\varepsilon} (t^{n+1}-t)} \hat{\mathcal{H}}_{p,l}(t^n) \, dt \right. \\
& \left. + \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\varepsilon} (t^{n+1}-t)} (t-t^n) \frac{d}{dt} \hat{\mathcal{H}}_{p,l}(t^n) \, dt \right) e^{ils} \\
& + \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\varepsilon} (t^{n+1}-t)} (t-t^n)^2 \left[ \int_0^1 (1-\theta) \frac{d^2}{dt^2} \hat{\mathcal{H}}_{p,l}(t^n + \theta(t-t^n)) \, d\theta \right] dt \right) e^{ils}.
\end{aligned}
 $$

Using $\sum_{l \in \mathbb{Z}} e^{ils} \hat{\mathcal{H}}_{p,l}(t^n) = \mathcal{H}_p(t^n, s)$ and a Taylor-Lagrange formula on $\frac{d}{dt} \hat{\mathcal{H}}_{p,l}(t^n)$, we get

$$ 
\begin{aligned}
\mathcal{Z}_p(t^{n+1}, s) = \mathcal{Z}_p(t^n, s - \frac{\Delta t}{\varepsilon}) &+ \int_{t^n}^{t^{n+1}} \mathcal{H}_p(t^n, s - \frac{t^{n+1}-t}{\varepsilon}) \, dt \\
&+ \int_{t^n}^{t^{n+1}} \frac{(t-t^n)}{\Delta t} \left( \mathcal{H}_p(t^n, s - \frac{t^{n+1}-t}{\varepsilon}) - \mathcal{H}_p(t^{n-1}, s - \frac{t^{n+1}-t}{\varepsilon}) \right) \, dt + r_p^n(s),
\end{aligned}
\qquad (55) $$
where $r_p^n(s)$ denotes the truncature error defined by

$$
\begin{align}
r_p^n(s) &= \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\epsilon}(t^{n+1}-t)} (t-t^n) \left[ \int_0^1 (1-\theta) \Delta t \frac{d^2}{dt^2} \hat{H}_{p,l}(t^n - \theta \Delta t) \, d\theta \right] dt \right) e^{ils} \nonumber \\
&\quad + \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\epsilon}(t^{n+1}-t)} (t-t^n)^2 \int_0^1 (1-\theta) \frac{d^2}{dt^2} \hat{H}_{p,l}(t^n + \theta(t-t^n)) \, d\theta \, dt \right) e^{ils} \tag{56} \\
&= \int_{t^n}^{t^{n+1}} \int_0^1 (1-\theta)(t-t^n) \Delta t \frac{d^2}{dt^2} \mathcal{H}_p(t^n - \theta \Delta t, s - \frac{t^{n+1}-t}{\epsilon}) \, d\theta \, dt \nonumber \\
&\quad + \int_{t^n}^{t^{n+1}} \int_0^1 (1-\theta)(t-t^n)^2 \frac{d^2}{dt^2} \mathcal{H}_p(t^n + \theta(t-t^n), s - \frac{t^{n+1}-t}{\epsilon}) \, d\theta \, dt. \nonumber
\end{align}
$$

Let us recall that the time derivative of $\hat{H}_{p,l}$ also involves time derivatives of $\mathcal{Z}_p$ and derivatives $\partial_z \hat{H}_{p,l}$. With the assumptions performed above, these quantities are all uniformly bounded. In the same way, the numerical scheme (37) on $\mathcal{Z}$ rewrites

$$
\begin{align}
\mathcal{Z}_p^{n+1}(s) &= \mathcal{Z}_p^n \left( s - \frac{\Delta t}{\epsilon} \right) + \sum_{l \in \mathbb{Z}} \left( \Delta t \varphi_0 \left( \frac{il}{\epsilon} \Delta t \right) \hat{H}_{p,l}^n + \Delta t^2 \varphi_1 \left( \frac{il}{\epsilon} \Delta t \right) \left[ \frac{\hat{H}_{p,l}^n - \hat{H}_{p,l}^{n-1}}{\Delta t} \right] \right) \nonumber \\
&= \mathcal{Z}_p^n \left( s - \frac{\Delta t}{\epsilon} \right) + \int_{t^n}^{t^{n+1}} \mathcal{H}_p^n \left( s - \frac{t^{n+1}-t}{\epsilon} \right) \, dt \tag{57} \\
&\quad + \int_{t^n}^{t^{n+1}} \frac{(t-t^n)}{\Delta t} \left[ \mathcal{H}_p^n \left( s - \frac{t^{n+1}-t}{\epsilon} \right) - \mathcal{H}^{n-1} \left( s - \frac{t^{n+1}-t}{\epsilon} \right) \right] \, dt, \nonumber
\end{align}
$$

where we recall that $\Delta t^{\alpha+1} \varphi_\alpha(\frac{il}{\epsilon}\Delta t) = \int_{t^n}^{t^{n+1}} e^{-\frac{il}{\epsilon}(t^{n+1}-t)} (t-t^n)^\alpha \, dt$ for $\alpha = 0, 1$.

Now, in view of deriving an induction on the error $e_{\mathcal{Z}_p}^n(s) = \mathcal{Z}_p^n(t^n, s) - \mathcal{Z}_p^n(s)$, we define, for $n \le N$

$$
\Delta \mathcal{H}_p^n(s) = \mathcal{H}_p \left( t^n, s - \frac{t^{n+1} - t}{\epsilon} \right) - \mathcal{H}_p^n \left( s - \frac{t^{n+1} - t}{\epsilon} \right).
$$

With this, considering the difference between (55) and (57), we obtain the following induction on $e_{\mathcal{Z}_p}^n(s)$

$$
\begin{equation}
\begin{split}
e_{\mathcal{Z}_p}^{n+1}(s) = e_{\mathcal{Z}_p}^n \left(s - \frac{\Delta t}{\epsilon}\right) &+ \int_{t^n}^{t^{n+1}} \Delta \mathcal{H}_p^n(s) \, dt \\
&+ \int_{t^n}^{t^{n+1}} \frac{(t-t^n)}{\Delta t} \left[ \Delta \mathcal{H}_p^n(s) - \Delta \mathcal{H}_p^{n-1}(s - \frac{\Delta t}{\epsilon}) \right] \, dt + r_p^n(s).
\end{split}
\tag{58}
\end{equation}
$$

In the following, we denote by $C$ all constants involved, even if it means taking the maximum of them. Constant $C$ does not depend on $\epsilon$ and $\tau$ because of (51). Under the assumptions on the nonlinear term $H$, we have

$$
|\Delta \mathcal{H}_p^n(s)| \leq C |e_{\mathcal{Z}_p}^n(s)| \leq C \|e_{\mathcal{Z}_p}^n\|_{L^\infty}.
$$

Moreover, the uniform estimates of $\partial_t^k \mathcal{Z}_p$ ($k = 0, 1, 2$) ensure that the local error can be estimated by

$$
|r_p^n(s)| \le \Delta t^3 \left\| \frac{d^2}{dt^2} \mathcal{H}_p \right\|_{L^\infty} \le C \Delta t^3. \quad (59)
$$
Setting $t_{\max} = N\Delta t$ and taking the maximum on $p$ leads to

$$ \|e_Z^{n+1}\|_{L^\infty} \le (1 + C\Delta t)\|e_Z^n\|_{L^\infty} + C\Delta t\|e_Z^{n-1}\|_{L^\infty} + C\Delta t^3. \quad (60) $$

Let us denote $U^n = (\|e_Z^n\|_{L^\infty}, \|e_Z^{n-1}\|_{L^\infty})^T$,

$$ A = \begin{pmatrix} 1 + C\Delta t & C\Delta t \\ 1 & 0 \end{pmatrix}, \quad b = \begin{pmatrix} C\Delta t^3 \\ 0 \end{pmatrix}. $$

We write (60) as $U^{n+1} \le AU^n + b$ to get recursively $U^n \le A^n U^0 + \sum_{k=0}^{n-1} A^k b$. Taking the norm, we get

$$ U^N \le (1 + 2C\Delta t)^N U^0 + \sum_{k=0}^{N-1} (1 + 2C\Delta t)^k b. $$

Once again, we denote by $C$ the constant involved. Since $U^0 = \max(\|e^0\|, \|e^{-1}\|) = \|e^{-1}\| = C\Delta t^2$ (iteration $n = -1$ computed with the order 1 scheme), we have finally

$$ U^N \le Ce^{CT}\Delta t^2 + e^{CT}\sum_{k=0}^{N-1} b \le Ce^{CT}\Delta t^2. $$

Let consider now the numerical scheme on the weights. First, using the same calculations (Taylor-Lagrange formula on $\hat{M}_{p,l}(t)$) as before, the two-scale solution $W_p$ of (52) satisfies

$$ \begin{aligned} W_p(t^{n+1}, s) = e^{-\frac{\Delta t}{\tau}} W_p(t^n, s - \frac{\Delta t}{\varepsilon}) &+ \frac{1}{\tau} \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-(\frac{il}{\varepsilon} + \frac{1}{\tau})(t^{n+1} - t)} \hat{M}_{p,l}(t^n) \, dt \right) e^{ils} \\ &+ \frac{1}{\tau} \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-(\frac{il}{\varepsilon} + \frac{1}{\tau})(t^{n+1} - t)} (t - t^n) \frac{\hat{M}_{p,l}(t^n) - \hat{M}_{p,l}(t^{n-1})}{\Delta t} \, dt \right) + q_p^n(s), \end{aligned} \quad (61) $$

where $q_p^n(s)$ denotes the truncature error

$$ \begin{aligned} q_p^n(s) = & \frac{1}{\tau} \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-(\frac{il}{\varepsilon} + \frac{1}{\tau})(t^{n+1} - t)} (t - t^n) \int_0^1 (1 - \theta) \Delta t \frac{d^2}{dt^2} \hat{M}_{p,l}(t^n - \theta \Delta t) \, d\theta \, dt \right) e^{ils} \\ & + \frac{1}{\tau} \sum_{l \in \mathbb{Z}} \left( \int_{t^n}^{t^{n+1}} e^{-(\frac{il}{\varepsilon} + \frac{1}{\tau})(t^{n+1} - t)} (t - t^n)^2 \int_0^1 (1 - \theta) \frac{d^2}{dt^2} \hat{M}_{p,l}(t^n + \theta(t - t^n)) \, d\theta \, dt \right) e^{ils}. \end{aligned} \quad (62) $$

Considering the difference between (61) and the scheme (37) on $W$ leads to an induction on the error

$$ \begin{aligned} e_{W_p}^{n+1}(s) = & e^{-\frac{\Delta t}{\tau}} e_{W_p}^n \left(s - \frac{\Delta t}{\varepsilon}\right) + \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-\frac{t^{n+1}-t}{\tau}} \Delta \mathcal{M}_p^n(s) \, dt \\ & + \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-\frac{t^{n+1}-t}{\tau}} \frac{t-t^n}{\Delta t} \, dt \left[ \Delta \mathcal{M}_p^n(s) - \Delta \mathcal{M}_p^{n-1}(s) \right] + q_p^n(s), \end{aligned} \quad (63) $$

where we introduced the following notation for $n \le N$

$$ \Delta \mathcal{M}_p^n(s) = \mathcal{M}_p \left( t^n, s - \frac{t^{n+1} - t}{\varepsilon} \right) - \mathcal{M}_p^n \left( s - \frac{t^{n+1} - t}{\varepsilon} \right). $$
Thanks to the definition of the $\varphi$ function (35) and the assumptions on $M$, the truncature error $q_p^n(s)$ given by (62) satisfies the following estimate

$$|q_p^n(s)| \le \left\| \frac{d^2 \mathcal{M}_p}{dt^2} \right\|_{L^\infty} \frac{1}{\tau} \left[ \int_{t^n}^{t^{n+1}} e^{-\frac{1}{\tau}(t^{n+1}-t)} (t-t^n) \Delta t \, dt + \int_{t^n}^{t^{n+1}} e^{-\frac{1}{\tau}(t^{n+1}-t)} (t-t^n)^2 \, dt \right] \\ \le \left\| \frac{d^2 \mathcal{M}_p}{dt^2} \right\|_{L^\infty} \frac{\Delta t^3}{\tau} \left[ \varphi_1 \left( \frac{\Delta t}{\tau} \right) + 2\varphi_2 \left( \frac{\Delta t}{\tau} \right) \right].$$

Consequently, under the assumptions on $M$ and the uniform estimates on $e_{Z_p}$, we have $|\Delta M_p^n(s)| \le C \|e_{Z_p}^n\|_{L^\infty} \le C\Delta t^2$. Then, the second term in the right hand side of (63) can be estimated by

$$\left| \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-\frac{t^{n+1}-t}{\tau}} \Delta M_p^n(s) \, dt \right| \le C \frac{\Delta t^3}{\tau} \varphi_0 \left( \frac{\Delta t}{\tau} \right),$$

whereas the third term satisfies

$$\left| \frac{1}{\tau} \int_{t^n}^{t^{n+1}} e^{-\frac{t^{n+1}-t}{\tau}} (t-t^n) \, dt \frac{\Delta M_p^n(s) - \Delta M_p^{n-1}(s)}{\Delta t} \right| \le 2C \frac{\Delta t^3}{\tau} \varphi_1 \left( \frac{\Delta t}{\tau} \right).$$

Thus, (63) becomes

$$\|e_W^{n+1}\|_{L^\infty} \le e^{-\frac{\Delta t}{\tau}} \|e_W^n\|_{L^\infty} + C \frac{\Delta t^3}{\tau} \left[ \varphi_0 \left( \frac{\Delta t}{\tau} \right) + 3\varphi_1 \left( \frac{\Delta t}{\tau} \right) + 2\varphi_2 \left( \frac{\Delta t}{\tau} \right) \right].$$

By induction, we have, with $N\Delta t = t_{\max}$,

$$\|e_W^N\|_{L^\infty} \le e^{-t_{\max}/\tau} \|e_W^0\|_{L^\infty} + C \frac{\Delta t^3 (1 - e^{-t_{\max}/\tau})}{\tau (1 - e^{-\Delta t/\tau})} (\varphi_0 + 3\varphi_1 + 2\varphi_2) \left(\frac{\Delta t}{\tau}\right), \\ \le e^{-t_{\max}/\tau} \|e_W^0\|_{L^\infty} + C(1 - e^{-t_{\max}/\tau}) \left[ \Delta t^2 + 2\tau^2 - 3\tau\Delta t + \frac{4\Delta t^2 - 2\tau\Delta t}{1 - e^{-\Delta t/\tau}} \right].$$

One can prove that $z^2 + 2 - 3z + \frac{4z^2 - 2z}{1 - e^{-z}} \le 5z^2$ for all $z > 0$ from which we deduce, using $\|e_W^0\|_{L^\infty} = 0$,

$$\|e_W^N\|_{L^\infty} \le e^{-t_{\max}/\tau} \|e_W^0\|_{L^\infty} + C(1 - e^{-t_{\max}/\tau})5\Delta t^2 \le C(1 - e^{-t_{\max}/\tau})\Delta t^2,$$

which concludes the proof. $\square$

**Remark 3.** When $M$ depends on $\omega$, in the truncature error (62) the term $\frac{d^2}{dt^2}M$ makes appear $\partial_t^2 W$. For a fixed $\tau$, time derivatives are uniformly bounded in $\varepsilon$, which enable to garantee the UA property in $\varepsilon$. However, when $\tau \to 0$ the time derivatives of $W$ are not uniformly bounded in $\tau$. In this case, we thus lose uniform accuracy in $\tau$ on weights (this property is still present on particles since it does not depend on $\tau$). For the AP property when $\tau \to 0$, we refer to Subsection 3.4.

# 5 Numerical results

This section is devoted to numerical illustrations of the properties of the schemes presented above. We first focus on a one-particle case and then on a PIC example where all the particles are coupled through the collision operator.

## 5.1 ODE framework

We consider the following multiscale ODE satisfied by $z = (x, v) \in \mathbb{R}^6$ and $\omega \in \mathbb{R}$

$$\begin{aligned} \dot{z}(t) &= \frac{1}{\varepsilon}Az(t) + h(t, z), \\ \dot{\omega}(t) &= \frac{1}{\tau}(M(\omega, z) - \omega(t)), \end{aligned} \qquad (64)$$
with $A$ and $h$ given by (3) and $E(t, x) = ((x_1 + x_3)\cos(t), x_1x_2\sin(t), -x_2^2e^{-t^2})$. Two examples for $\mathcal{M}$ will be investigated. In a first example, $\mathcal{M}(\omega, z) = 1 + |z|^2 + e^{-\omega^2} (\mathcal{M}$ depends on $\omega$) whereas in a second example $\mathcal{M}(z) = 1 + |z|^2$ ($\mathcal{M}$ does not depend on $\omega$). The initial condition is $z(0) = (1, 1, 0, 1/2, 1/2, 3/2), \omega(0) = 1$ and the final time of the simulation is $t_{\max} = 1$.

In order to illustrate the properties of the schemes, we plot the error between the numerical and the exact solutions. As we do not know the exact one, we replace it by a reference solution obtained with a small time step $\Delta t = 2^{-14}$. The error considered is the sum of the $L^\infty$-error of the solution components. In all tests, we consider the second order scheme (37)-(38)-(39) with the initial condition (41) and $N_s = 16$ are used to discretize the $s$ variable.

In Figures 1, we consider the case $\mathcal{M}(\omega, z)$. We plot the error as a function of the time step $\Delta t = 2^{-k}, k = 3, \dots, 13$ for different values of $\varepsilon = 10^{-k}(k = 0, \dots, 6)$ with $\tau = 1$ (left figure) and of $\tau = 10^{-k}(k = 0, \dots, 4)$ with $\varepsilon = 1$ (right figure). We observe the uniform accuracy with respect to $\varepsilon$ (left figure) since for all the values of $\varepsilon$, the error remains the same and the second order is recovered, as expected. However, it is not the case when $\tau$ varies (right figure) since the order is observed only when $\tau > \Delta t$. In Figure 2, we consider the case $\mathcal{M}(z)$ for which Theorem 1 ensures uniform accuracy in both $\varepsilon$ and $\tau$ parameters. As expected, we obtain the uniform accuracy property in $\tau$ (the uniform accuracy in $\varepsilon$ is ensured by the two-scale strategy).

<figure><img src="image_2.png" /><figcaption>Figure 1: Case $\mathcal{M}(w, z)$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Case $\mathcal{M}(z): L^\infty$-Error for $\varepsilon = 1$ and different $\tau$.</figcaption></figure>

In Figures 3, we plot the time history of the Maxwellian weights $\mathcal{M}(z(t))$ together with $w(t)$ for different values
of $\epsilon$ and $\tau$, with $\Delta t = 2.5 \times 10^{-4}$. These two figures illustrate the asymptotic properties (gyrokinetic limit $\epsilon \to 0$ and fluid limit $\tau \to 0$) of (64). Indeed, when $\epsilon$ is fixed, we observe the convergence when $\tau \to 0$ of $w(t)$ to $M(z(t))$ which oscillates in time (see Figure 3-(a)). However, for $\tau > 0$ fixed, and $\epsilon \to 0$, we observe the convergence of $w(t)$ towards the average of $M(z(t))$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Time history of the $M(z(t))$ and of $\omega(t)$ for different values of $\epsilon$ and $\tau$.</figcaption></figure>

Finally, in Figure 4, we illustrate the Uniform Accuracy property in both parameters $\epsilon$ and $\tau$. We take $\epsilon = \tau = 0.01$ and plot the time history of $\omega$ for two different values of the time step: $\Delta t = 2.5 \times 10^{-4}$ and $\Delta t = 0.05$. We can observe that even if $\Delta t$ is bigger, the numerical scheme is able to capture the strong relaxation (at the very beginning of the simulation) and the highly oscillatory behavior, without resolving these two stiffnesses.

<figure><img src="image_5.png" /><figcaption>Figure 4: Weights trajectories for $\Delta t$ bigger than $\tau$.</figcaption></figure>

## 5.2 PDE framework

In this part, we consider the PIC discretization of (2)-(3) with a BGK operator. The goal is to observe the relaxation to the gyromaxwellian equilibria. Since the problem is inherently six-dimensional and as such, is very costly from a numerical point of view, we simplify the problem by considering $h = 0$ in (2)-(3) and remove the dependence with respect to $x_3$. As a consequence, the filtered model satisfied by $F^{\epsilon, \tau}(t, Z)$ becomes

$$ \partial_t F^{\epsilon, \tau} = \frac{1}{\tau} Q_{\text{filt}}[F^{\epsilon, \tau}](t, t/\epsilon, Z) = \frac{1}{\tau} (M_{\text{filt}}[F^{\epsilon, \tau}] - F^{\epsilon, \tau}), $$
where $M_{\text{filt}}$ is defined in (12). It follows that the gyromoments, and thus the equilibria, are independent of time. The goal is to observe the relaxation towards the gyromaxwellian in the gyrofluid limit ($\epsilon \to 0$ then $\tau \to 0$). Hence, from an out-of-equilibria initial data $F_0$, we compute its associated equilibria $\mathcal{G}[F_0]$, and observe the relaxation of $F^{\epsilon, \tau}$ towards $\mathcal{G}[F_0]$. For the initial condition $F_0$, we choose it as the sum of two different gyromaxwellians of same mass characterized respectively by the following macroscopic quantities $U_a = (n, U_a, \bar{Y}_a, K_a, G_a)$ and $U_b = (n, U_b, \bar{Y}_b, K_b, G_b)$. After some computations, we find that the gyromaxwellian $\mathcal{G}[F_0]$ associated to $F_0$ is characterized by

$$ U_{ab} = \left( 2n, \frac{U_a+U_b}{2}, \frac{\bar{Y}_a+\bar{Y}_b}{2}, \frac{K_a+K_b}{2} + \frac{|U_a-U_b|^2}{8}, \frac{G_a+G_b}{2} + \frac{|\bar{Y}_a-\bar{Y}_b|^2}{8} - \frac{|\bar{U}_a-\bar{U}_b|^2}{8} \right). $$

For the first test, we consider the scheme B (36a)-(47)-(48) with the following numerical parameters: $N_p = 1.6 \times 10^7$ particles, $N_s = 8$ for the s direction and $\Delta t = 2 \times 10^{-4}$ whereas $\epsilon = 10^{-8}$ and $\tau = 10^{-4}$.

In Figure 5, we illustrate the relaxation from a sum of two gyromaxwellians of gyromoments $U_a = (0.5, 0, 0, 0, -1, 0, 0.3, 0)$ and $U_b = (0.5, 0, 0, 0, 1, 0, 0.3, 0)$ to a gyromaxwellian $\mathcal{G}[F_0]$ of gyromoments $U_{ab} = (1, 0, 0, 0, 0, 0, 0.3, 0.5)$. More precisely, from the particle approximation of $F^{\epsilon, \tau}(t, Z)$, we plot a reconstruction of $F^{\epsilon, \tau}(t = 0, X_1, X_2 = 0, V_1 = 0, V_2 = 0, V_3 = 0)$ and $F^{\epsilon, \tau}(t = t_{\text{max}}, X_1, X_2 = 0, V_1 = 0, V_2 = 0, V_3 = 0)$ (let us recall that the $X_3$-dimension has been removed). These $X_1$-slices are plotted in Figure 5 together with the analytical expression of $F_0$ and $\mathcal{G}[F_0]$ slices. Let us remark that computing a one-dimensional slice of a particles ensemble requires some interpolation operators. Here, before taking slice we have regularized the particle density in space and velocity variables with the same B-spline as before. Even if the method has some noise, we clearly observe the relaxation to the correct limit represented by the gyromaxwellian.

<figure><img src="image_6.png" /><figcaption>Figure 5: Particle density relaxing to a gyromaxwellian.</figcaption></figure>

We are now interested in more quantitative results by considering the time evolution of the error between the particles weights and the gyromaxwellian weights. To do so, we plot at each time iteration the three following errors: $e_1^n = \sum_p |\omega_p^n - \mathcal{G}_p^{ex}|$, $e_2^n = \sum_p |\omega_p^n - \mathcal{G}_p^n|$ and $e_3^n = \sum_p |\mathcal{G}_p^n - \mathcal{G}_p^{ex}|$, where $\mathcal{G}_p^{ex}$ denotes the weights computed from analytical gyromaxwellian $\mathcal{G}[F_0]$ (which are constant in time) and $\mathcal{G}_p^n$ are numerical gyromaxwellian weights computed as in part 3.5. These quantities are plotted for the two schemes introduced above: the scheme A (36a)-(44) and the scheme B (36a)-(47)-(48). The parameters taken are $N_p = 12288000$, $\epsilon = 10^{-8}$, $\Delta t = 10^{-3}$, $N_s = 4$, $\Delta x \approx 0.37$ and we look at different values of $\tau$.

In Figures 6 and 7, we can see that $e_3^n$ increases in time whereas it should remain small. This is due to the fact that gyromoments are not conserved exactly by our scheme which deteriorates the computation of the gyromaxwellian $\mathcal{G}_p^n$. A first way to improve this conservation is to increase the number of particles (we checked that the conservations are improved by adding particles). A second option would be to use a projection technique as discussed in [16, 21] to enforce the conservation of the gyromoments at the discrete level.

Let us now focus on the other errors $e_1^n$ and $e_2^n$ for the two schemes. In Figure 6, the error $e_2^n$ computed with the scheme A does not present the good behavior since the convergence towards the gyromaxwellian is too slow. However, we can observe in Figure 7 that the scheme B is able to capture the right limit in one iteration (since $\Delta t > \tau$) thanks to its AP property. The error $e_1^n$ is the error $e_2^n$ with the effect of non conservation of gyromoments which are deteriorated with time.
<figure><img src="image_7.png" /><figcaption>Figure 6: Evolution of the three errors $e_i^n$, $i = 1, 2, 3$ for scheme A for two different $\tau$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Evolution of the three errors $e_i^n$, $i = 1, 2, 3$ for scheme B for two different $\tau$.</figcaption></figure>

# 6 Conclusion and opening

In this paper, we have constructed and analysed multiscale numerical schemes for the collisional Vlasov equation in the finite Larmor radius approximation regime. This model, proposed in [4], involves two scaling parameters: $\epsilon$ drives the oscillation due to the presence of the strong external magnetic field and $\tau$ drives the hydrodynamic limit. We also made an analysis of the gyrofluid limit model.

Further developments can be performed. First, regarding the computational cost (due to high dimensionality of the problem), it would be interesting to investigate a micro-macro approach where the one-dimensional gyrofluid model would be coupled to a particle discretization of a collisional Vlasov equation, in the spirit of [11, 27]. Considering the coupling with the self-consistent electromagnetic field governed by Maxwell equation could be also considered.
# Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

This work has been also supported by ANR project MUFFIN (ANR-19-CE46-0004) and by Centre Henri Lebesgue, programme ANR-11-LABX-0020-0.

## A Gyrofluid model's properties

In this section, we study the gyrofluid model derived from the gyrokinetic model (8) with the BGK operator when $\tau \to 0$. The main interest of this gyrofluid model lies in the fact that the macroscopic unknown only depend on time and on one spatial variable $X_3$. For the sake of simplicity, we assume here that the electric field is given and only depends on time: $E = E(t)$ but the following calculations can be generalized to a given space dependent electric field. Considering the gyromoments of (8) and replacing $F^\tau$ by the gyromaxwellian (15) enables to get the following gyrofluid model satisfied by $U(t, X_3) = (n, n\bar{U}, nU_3, n\bar{Y}, n(K + |U|^2/2), n(G + |\bar{Y}|^2/2 - |\bar{U}|^2/2))(t, X_3) \in \mathbb{R}^8$

$$ \partial_t n + \partial_{X_3} n U_3 = 0, \qquad (65a) $$

$$ \partial_t n \bar{U} + \partial_{X_3} n U_3 \bar{U} = 0, \qquad (65b) $$

$$ \partial_t n U_3 + \partial_{X_3} (n U_3^2 + n \theta) = n E_3, \qquad (65c) $$

$$ \partial_t n \bar{Y} + \partial_{X_3} n U_3 \bar{Y} = n \perp \bar{E}, \qquad (65d) $$

$$ \partial_t \left[ n \left( K + \frac{|U|^2}{2} \right) \right] + \partial_{X_3} \left[ n U_3 \left( K + \frac{|U|^2}{2} + \theta \right) \right] = n U_3 E_3, \qquad (65e) $$

$$ \partial_t \left[ n \left( G + \frac{|\bar{Y}|^2}{2} - \frac{|\bar{U}|^2}{2} \right) \right] + \partial_{X_3} \left[ n U_3 \left( G + \frac{|\bar{Y}|^2}{2} - \frac{|\bar{U}|^2}{2} \right) \right] = n \bar{Y} \cdot \perp \bar{E}. \qquad (65f) $$

The gyrofluid model has also an entropic inequality, derived directly from H-theorem's inequality

$$ \partial_t n \ln \left( \frac{n}{\frac{\mu^2 \theta^{3/2}}{\mu - \theta}} \right) + \partial_{X_3} n U_3 \ln \left( \frac{n}{\frac{\mu^2 \theta^{3/2}}{\mu - \theta}} \right) \leq 0, \qquad (66) $$

$\mu$ and $\theta$ being defined by (18). A space dependent electric field would have given a different right hand side in (65): the change of variable (4) and the averaging have to be applied to $E$ to get an average electric field as a source term. One can observe that the gyrofluid model (65) shares similarities with compressible Euler equations but also has its own specificities. Indeed, in addition to the equation of continuity, there are five equations on mean velocity (with $U_3$ as the main velocity), and two equations (65e) and (65f) on the kinetic and gyration energies, the sum of which gives the equation on the total energy. The kinetic (resp. gyration) energy is composed of an internal kinetic (resp. gyration) energy $K$ (resp. $G$) and a kinetic (resp. gyration) energy $|U|^2/2$ (resp. $(|\bar{Y}|^2 - |\bar{U}|^2)/2$). Moreover, as $\theta$ is a function of the two internal energies $K$ and $G$, $n\theta$ can be interpreted as an equivalent of the pressure in the compressible Euler equations. The right hand side terms correspond to external electric field applied on the system. The gyrofluid model (65) is presented in a conservative form with a source term. In order to study algebra and wave-structure of the model, we rewrite it in a quasi-linear form expressed with primitive variables $W(t, X_3) = (n, U, \bar{Y}, K, G)(t, X_3) \in \mathbb{R}^8$. This is done in the following Proposition.

**Proposition 2.** The gyrofluid model (65) can be rewritten in a quasilinear form

$$ \partial_t W + A(W) \partial_{X_3} W = B, \qquad (67) $$
with $\mathbf{W}(t, X_3) = (n, U, \bar{Y}, K, G)(t, X_3)$ and with

$$ \mathbf{A}(\mathbf{W}) = \begin{pmatrix} U_3 & 0 & 0 & n & 0 & 0 & 0 & 0 \\ 0 & U_3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & U_3 & 0 & 0 & 0 & 0 & 0 \\ \frac{\theta}{n} & 0 & 0 & U_3 & 0 & 0 & \theta_K & \theta_G \\ 0 & 0 & 0 & 0 & U_3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & U_3 & 0 & 0 \\ 0 & 0 & 0 & \theta & 0 & 0 & U_3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & U_3 \end{pmatrix}, \qquad \mathbf{B} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ E_3(t) \\ E_2(t) \\ -E_1(t) \\ 0 \\ 0 \end{pmatrix}. \qquad (68) $$

Here $\theta_K$ and $\theta_G$ are such that $\partial_{X_3}\theta = \theta_K\partial_{X_3}K + \theta_G\partial_{X_3}G$, $\theta$ being defined by (18). The gyrofluid model is hyperbolic: $\mathbf{A}(\mathbf{W})$ is diagonalizable and presents the eigenvalue $\lambda_0 = U_3$ of multiplicity 6 and two eigenvalues $\lambda_- = U_3 - c$ and $\lambda_+ = U_3 + c$ with $c = \sqrt{\theta(1 + \theta_K)}$ of multiplicity 1. The associated eigenvectors for $\lambda_0$ are

$$ \begin{aligned} r_1 &= (0, 1, 0, 0, 0, 0, 0, 0), & r_2 &= (0, 0, 1, 0, 0, 0, 0, 0), & r_3 &= (0, 0, 0, 0, 1, 0, 0, 0), \\ r_4 &= (0, 0, 0, 0, 0, 1, 0, 0), & r_5 &= (0, 0, 0, 0, 0, 0, \theta_G, -\theta_K), & r_6 &= (-n\theta_K, 0, 0, 0, 0, 0, \theta, 0). \end{aligned} \qquad (69) $$

The other eigenvectors are given by

$$ r_{-} = (n, 0, 0, -c, 0, 0, \theta, 0) \text{ for } \lambda_{-}, \quad r_{+} = (n, 0, 0, +c, 0, 0, \theta, 0) \text{ for } \lambda_{+}. $$

Remark that similarly to the classical Euler equations, the gyrofluid model presents two genuinely nonlinear fields associated to $\lambda_{\pm} = U_3 \pm c$, as well as a linearly degenerated field associated to $\lambda_0 = U_3$. To prove Proposition (2), we need the following Lemma.

**Lemma 1.** We have the relations below

$$ \partial_t n \frac{|\bar{U}|^2}{2} + \partial_{X_3} n U_3 \frac{|\bar{U}|^2}{2} = 0, \qquad (70) $$

$$ \partial_t n \frac{|\bar{Y}|^2}{2} + \partial_{X_3} n U_3 \frac{|\bar{Y}|^2}{2} = n \bar{Y} \cdot \perp \bar{E}, \qquad (71) $$

$$ \partial_t n \frac{|U|^2}{2} + \partial_{X_3} n U_3 \frac{|U|^2}{2} = n U_3 E_3 - U_3 \partial_{X_3} n \theta. \qquad (72) $$

*Proof.* Starting from (65b) and use the continuity equation (65a) enables to get $\partial_t \bar{U} + U_3 \partial_{X_3} \bar{U} = 0$. Consequently, developing (70) leads to

$$ \partial_t n \frac{|\bar{U}|^2}{2} + \partial_{X_3} n U_3 \frac{|\bar{U}|^2}{2} = \frac{|\bar{U}|^2}{2} (\partial_t n + \partial_{X_3} n U_3) + n(\bar{U} \cdot (\partial_t \bar{U} + U_3 \partial_t \bar{U})) = 0. $$

Relations (71) and (72) are obtained by similar calculations.

We prove now Proposition 2.

*Proof.* Using Lemma 1 together with some calculations, we will show that (65) can be written into (67) for which the eigenstructure can be obtained easily. First, (65a) and (65b) give

$$ \partial_t n + U_3 \partial_{X_3} n + n \partial_{X_3} U_3 = 0 \quad \text{and} \quad \partial_t \bar{U} + U_3 \partial_{X_3} \bar{U} = 0, $$

whereas (65c) gives

$$ \partial_t U_3 + U_3 \partial_{X_3} U_3 + \partial_{X_3} \theta + \frac{\theta}{n} \partial_{X_3} n = E_3, $$

where $\theta$ is a function of K and G which are independent of the other variables (see (18) and (19)). Then, there exist $\theta_K(K, G)$ and $\theta_G(K, G)$ such that $\partial_{X_3}\theta = \theta_K\partial_{X_3}K + \theta_G\partial_{X_3}G$. Consequently, (65c) becomes

$$ \partial_t U_3 + U_3 \partial_{X_3} U_3 + \theta_K \partial_{X_3} K + \theta_G \partial_{X_3} G + \frac{\theta}{n} \partial_{X_3} n = E_3. \qquad (73) $$
From (65d), we obtain

$$ \partial_t \bar{Y} + U_3 \partial_{X_3} \bar{Y} = \perp \bar{E}. $$

We use (72) of Lemma 1 together with (65e) to get

$$ \partial_t n K + \partial_{X_3} n U_3 K + n \theta \partial_{X_3} U_3 = 0, $$

or again

$$ \partial_t K + U_3 \partial_{X_3} K + \theta \partial_{X_3} U_3 = 0. \qquad (74) $$

We use (70) and (71) of Lemma 1 to recast (65f) into

$$ \partial_t G + U_3 \partial_{X_3} G = 0. $$

We then deduce that (65) can be written in the quasilinear form (67). Moreover, we can check that $\theta_K > 0$ which implies that $c = \sqrt{\theta(1 + \theta_K)}$ is well defined. Because of the shape of A, we see that the characteristic polynomial is reduced to the computation of the determinant of a 3x3 matrix

$$ \det(A - \lambda I d) = (U_3 - \lambda)^5 \det \begin{pmatrix} U_3 - \lambda & n & 0 \\ \theta/n & U_3 - \lambda & \theta_K \\ 0 & \theta & U_3 - \lambda \end{pmatrix}. $$

It follows that $\det(A - \lambda I d) = (U_3 - \lambda)^6(\lambda^2 - 2U_3\lambda + U_3^2 - c^2) = (U_3 - \lambda)^6(\lambda - (U_3 - c))(\lambda - (U_3 + c))$. In the same way, we see that $r_1, r_2, r_3$ and $r_4$ are eigen vectors associated to $U_3$. Then, it remains to find other eigen vectors on the 4x4 matrix

$$ \begin{pmatrix} U_3 & n & 0 & 0 \\ \theta/n & U_3 & \theta_K & \theta_G \\ 0 & \theta & U_3 & 0 \\ 0 & 0 & 0 & U_3 \end{pmatrix}. $$

This concludes the proof. □

## References

[1] C. K. Birdsall and A. B. Langdon. *Plasma physics via computer simulation*. Taylor and Francis, New York, 2005.

[2] M. Bostan. High magnetic field equilibria for the Fokker-Planck-Landau equation. *Annales de l'Institut Henri Poincaré (C) Non Linear Analysis*, page 10.1016/j.anihpc.2015.01.008, 2015.

[3] M. Bostan and A. Finot. The Effective Vlasov-Poisson System for the Finite Larmor Radius Regime. *Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal*, 14(4):1238-1275, 2016.

[4] M. Bostan and A. Finot. The finite Larmor radius regime: collisional setting and fluid models. *Communications in Contemporary Mathematics*, 2019.

[5] P. Chartier, N. Crouseilles, M. Lemou, and F. Méhats. Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations. *Numerische Mathematik*, 129(2):211-250, 2015.

[6] P. Chartier, N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field. *Mathematics of Computation*, 88(320):2697-2736, 2019.

[7] P. Chartier, N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate methods for three dimensional Vlasov equations under strong magnetic field with varying direction, 2019.

[8] P. Chartier, N. Crouseilles, and X. Zhao. Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime. *Journal of Computational Physics*, 375:619-640, Dec. 2018.

[9] P. Chartier, M. Lemou, F. Méhats, and X. Zhao. Derivative-free high-order uniformly accurate schemes for highly-oscillatory systems. *IMA Journal of Numerical Analysis*, 2022.
[10] F. Coron and B. Perthame. Numerical Passage from Kinetic to Fluid Equations. *SIAM Journal on Numerical Analysis*, 28(1):26-42, 1991.

[11] A. Crestetto, N. Crouseilles, and M. Lemou. Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles. *Kinetic and Related Models*, 5(4):787-816, 2012.

[12] A. Crestetto, N. Crouseilles, and M. Lemou. A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling. *Communications in Mathematical Sciences*, 16(4):887-911, 2018.

[13] N. Crouseilles, G. Dimarco, and M.-H. Vignal. Multiscale schemes for the BGK-Vlasov-Poisson system in the quasi-neutral and fluid limits. Stability analysis and first order schemes. *Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal*, 14(1):65-95, 2016.

[14] N. Crouseilles, M. Lemou, and F. Méhats. Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations. *Journal of Computational Physics*, 248:287-308, 2013.

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

[16] G. Dimarco and R. Loubere. Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation. *Journal of Computational Physics*, 255:680-698, Dec 2013.

[17] G. Dimarco and L. Pareschi. Exponential Runge-Kutta Methods for Stiff Kinetic Equations. *SIAM Journal on Numerical Analysis*, 49(5):p. 2057-2077, 2011.

[18] F. Filbet and S. Jin. A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources. *Journal of Computational Physics*, 229(20):7625-7648, 2010.

[19] F. Filbet and S. Jin. An Asymptotic Preserving Scheme for the ES-BGK model. *Journal Scientific Computing*, 46(2):204-224, 2011.

[20] E. Frénod and E. Sonnendrücker. The Finite Larmor Radius Approximation. *SIAM Journal on Mathematical Analysis*, 32(6):1227-1247, 2001.

[21] I. M. Gamba and S. H. Tharkabhushanam. Spectral - Lagrangian methods for Collisional Models of Non - Equilibrium Statistical States, 2008.

[22] M. Hochbruck and A. Ostermann. Exponential integrators. *Acta Numer.*, 19:209-286, may 2010.

[23] D. Issautier. Convergence of a weighted particle method for solving the boltzmann (bgk) equation. *SIAM journal on numerical analysis*, 33(6):2099-2119, 1996.

[24] S. Jin. Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations. *SIAM J. Sci. Comput.*, 21(2):441-454, Sept. 1999.

[25] S. Jin. Asymptotic-Preserving Schemes for Multiscale Physical Problems. *Acta Numerica*, 2021.

[26] S. Jin and B. Yan. A Successive Penalty-Based Asymptotic-Preserving Scheme for Kinetic Equations. *SIAM Journal on Scientific Computing*, 35(1):A150-A172, 2013.

[27] M. Lemou and L. Mieussens. A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit. *SIAM Journal on Scientific Computing*, 31(1):334-368, 2008.

[28] Q. Li and J. Lu. An Asymptotic Preserving Method for Transport Equations with Oscillatory Scattering Coefficients. *Multiscale Modeling & Simulation*, 15(4):1694-1718, 2017.

[29] Y. Mocquard, P. Navaro, and N. Crouseilles. HOODESolver.jl: A Julia package for highly oscillatory problems. *J. Open Source Softw.*, 6:3077, 2021.

[30] W. Ren, H. Liu, and S. Jin. An asymptotic-preserving Monte Carlo method for the Boltzmann equation. *Journal of Computational Physics*, 276:380-404, 2014.
[31] E. Sonnendrücker. *Numerical Methods for the Vlasov-Maxwell equations*. Lecture notes, 2015.

[32] B. Wang and X. Zhao. Geometric two-scale integrators for highly oscillatory system: uniform accuracy and near conservations, 2021.

[33] B. Zhang, H. Liu, and S. Jin. An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation. *Journal of Computational Physics*, 305:575–588, 2016.