domain: math math.math-na phys

text:
Anaïs Crestetto

Nicolas Crouseilles

Damien Prel

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

Keywords: Vlasov equation, Finite Larmor radius, Multiscale numerical scheme MSC codes. 65M25, 65M75, 35Q83

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.

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 [START_REF] Bostan | The Effective Vlasov-Poisson System for the Finite Larmor Radius Regime[END_REF][START_REF] Frénod | The Finite Larmor Radius Approximation[END_REF]). 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 [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF] and [START_REF] Bostan | High magnetic field equilibria for the Fokker-Planck-Landau equation[END_REF]). Indeed, under some quite general assumptions, the authors proved in [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF] 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)

∂ t f + 1 ε (v 1 ∂ x1 f + v 2 ∂ x2 f ) + v 3 ∂ x3 f + E • ∇ v f + 1 ε (v 2 ∂ v1 f -v 1 ∂ v2 f ) = 1 τ Q(f ), (1)
where (t, x, v) ∈ R + × R 3 × R 3 (x = (x 1 , x 2 , x 3 ) and v = (v 1 , v 2 , v 3 )), E = (E 1 , E 2 , E 3 ) : R 3 → 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 [START_REF] Birdsall | Plasma physics via computer simulation[END_REF], which can run for arbitrary values of ε (scaled cyclotronic period) and τ (Knudsen number). Let us remark that the limit ε → 0 (τ fixed) leads to highly-oscillatory (in time) behavior of f whereas the limit τ → 0 (ε 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 [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF], using two-scale strategies developed 1 in [5-7, 9, 14, 15]. On the other side, Asymptotic Preserving numerical schemes efficient in the hydrodynamic limit have been widely studied [START_REF] Coron | Numerical Passage from Kinetic to Fluid Equations[END_REF][START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Dimarco | Exponential Runge-Kutta Methods for Stiff Kinetic Equations[END_REF][START_REF] Filbet | An Asymptotic Preserving Scheme for the ES-BGK model[END_REF][START_REF] Jin | Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations[END_REF][START_REF] Jin | Asymptotic-Preserving Schemes for Multiscale Physical Problems[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. However, when two limits are combined, the literature is less abundant (see [START_REF] Crouseilles | Multiscale schemes for the BGK-Vlasov-Poisson system in the quasi-neutral and fluid limits. Stability analysis and first order schemes[END_REF][START_REF] Li | An Asymptotic Preserving Method for Transport Equations with Oscillatory Scattering Coefficients[END_REF]).

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 [START_REF] Birdsall | Plasma physics via computer simulation[END_REF]) assumes the distribution function is approximated by a sum of macro-particles whose positions x(t), velocities v(t), and weights ω(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 [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Chartier | Derivative-free high-order uniformly accurate schemes for highly-oscillatory systems[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF]. 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 [START_REF] Birdsall | Plasma physics via computer simulation[END_REF], the particles weights ω(t) are solution to an equation of the form ω(t) = 1 τ (M (t/ε, ω(t)) -ω(t)) with M (θ, •) a nonlinear function of ω. 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 [START_REF] Dimarco | Exponential Runge-Kutta Methods for Stiff Kinetic Equations[END_REF][START_REF] Jin | A Successive Penalty-Based Asymptotic-Preserving Scheme for Kinetic Equations[END_REF], it is possible to design a scheme which enjoys the uniformly accuracy property with respect to ε (oscillations) and the asymptotic preserving property with respect to τ (relaxation). Moreover, in a reduced case where the filtered Maxwellian does not depend on the weigths, it is possible to prove that the numerical scheme enjoys the uniform accuracy in both ε and τ .

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.

Presentation of the models

In this section, the different asymptotic models of (1) are recalled from [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF]. We first recall the highly oscillatory limit ε → 0 (τ fixed): after filtering out the main oscillation, a collisional gyrokinetic model is obtained by a standard averaging procedure. The main point highlighted in [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF] lies in the fact that the averaged collision operator ⟨Q⟩ 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.

Gyrokinetic model

To formally derive the collisional gyrokinetic model from the limit ε → 0 of (1), we make explicit the ε and τ dependence in the distribution function unknown by denoting it f ε,τ (t, z) with z = (x, v), so that (1) is recast as

∂ t f ε,τ + 1 ε Az • ∇ z f ε,τ + h(t, z) • ∇ z f ε,τ = 1 τ Q(f ε,τ ), (2)
with

A =         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         , h(t, z) =         0 0 v 3 E 1 (t, x) E 2 (t, x) E 3 (t, x)         . ( 3
)
On the one side, the perpendicular (to the magnetic field) dynamics 1 ε 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 -t ε A z, (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 Z associated to Az • ∇ z by the relation Z(s, Z) = e sA Z. The flow is given by

Z = (X, V), X(s, Z) = (X + (I 2 -R(-s)) ⊥ V , X 3 ), V(s, Z) = (R(-s)V , V 3 ), (5)
where

X = (X 1 , X 2 ), V = (V 1 , V 2 ), ⊥ V = (V 2 , -V 1 )
, I 2 the 2x2 identity matrix and R(θ) stands for the rotation of angle θ. Due to the spectral properties of A, one can prove that the flow s → Z(s, Z) is 2π-periodic, which will have important consequences from a numerical point of view. Let now define F ε,τ the filtered distribution function F ε,τ (t, Z) = f ε,τ (t, z) (with z = e t/εA Z from ( 4)) which satisfies the filtered version of (2)

∂ t F ε,τ + H(t, t/ε, Z) • ∇ Z F ε,τ = 1 τ Q filt [F ε,τ ](t, t/ε, Z), (6)
where

Q filt [F ε,τ ](t, s, Z) = Q[F ε,τ (t, e -sA •)](e sA Z), H(t, s, Z) = e -sA h(t, e sA Z). (7)
As one can see, the definition of Q 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 filt . The filtered model ( 6) highlights three time scales. First, the slow time scale corresponding to t → 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 θ → H(t, θ/ε, Z); lastly, the presence of the collision term induces a collisional time scale t/τ . The averaged model obtained when ε goes to 0 is derived by averaging with respect to the fast time variable t/ε. Following [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF], the collisional gyrokinetic model satisfied by F τ (which is the limit as ε → 0 of F ε,τ ) writes

∂ t F τ + ⟨H⟩ • ∇ Z F τ = 1 τ ⟨Q⟩ [F τ ], (8)
with

⟨H⟩ (t, Z) = 1 2π 2π 0 H(t, s, Z) ds, ⟨Q⟩ [F τ ](t, Z) = 1 2π 2π 0 Q filt [F τ ](t, s, Z) ds. (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 [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF]) of the equilibrium and collision invariants of ⟨Q⟩ is necessary since, as τ → 0, F τ converges to an equilibrium of ⟨Q⟩.

The average operator ⟨Q BGK ⟩

In view of numerical simuations, we will use the standard BGK operator which is the easiest operator satisfying the properties required by Theorem 4.1 of [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF]. Indeed, the analysis presented in [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF] 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 BGK (f ) = M[f ] -f, (10)
where the maxwellian M[f ] associated to f is defined by

M[f ](x, v) = n(x) (2πT (x)) 3/2 exp - |v -u(x)| 2 2T (x) , (11)
where the moments n, u and T are defined by

n(x) = R 3 f (x, v) dv, nu(x) = R 3 f (x, v)v dv, n |u| 2 2 + 3 2 T (x) = R 3 f (x, v) |v| 2 2 dv.
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))

M filt [F ](s, Z) = M[F (e -sA •)](e sA Z), (12)
and the averaged Maxwellian by

⟨M⟩ [F ](Z) = 1 2π 2π 0 M filt [F ](s, Z) ds. ( 13
)
We give now the expression of the average operator

⟨Q BGK ⟩ ⟨Q BGK ⟩ [F ](Z) = ⟨M⟩ [F ](Z) -F (Z). ( 14
)
Thanks to Proposition 4.5 of [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF], the equilibria of this averaged BGK operator are the so-called gyromaxwellians defined by (recalling X = (X 1 , X 2 ) and

V = (V 1 , V 2 )) G[F ](X, X 3 , V ) = nM 2 µθ µ-θ (V -U )M 1 θ (V 3 -U 3 )M 2 µ (X + ⊥ V -Y ), (15)
with

M d T (v) = (2πT ) -d/2 exp -|v| 2 2T
where n(X

3 ) ∈ R + , U (X 3 ) ∈ R 3 , Y (X 3 ) ∈ R 2 , θ(X 3 ) ∈ R + and µ(X 3 ) ∈ R +
are defined by integrating F against the 8 collisional invariants of ⟨Q BGK ⟩ (whereas

Q BGK has 5 collisional invariants) n(X 3 ) = R 2 R 3 F (X, V ) dV dX, (16a)
nU (X 3 ) = R 2 R 3 F (X, V )V dV dX, nY (X 3 ) = R 2 R 3 F (X, V )(X + ⊥ V ) dV dX, ( 16b
) n K + |U | 2 2 (X 3 ) = R 2 R 3 F (X, V ) |V | 2 2 dV dX, ( 16c
)
n G+ Y 2 2 - U 2 2 (X 3 ) = R 2 R 3 F (X, V ) X + ⊥ V 2 -V 2 2 dV dX. ( 16d
)
These quantities, which only depend on X 3 , are called gyromoments: n is the density of particules, U = (U 1 , U 2 , U 3 ) is the average velocity, 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 θ and µ then verify µ > θ > 0 and satisfy the following relations

K = µθ µ -θ + θ 2 , G = µ - µθ µ -θ . ( 17
)
We can inverse the system to express θ and µ as a function of K and G to get

θ = 2(K + G)/(1 + 2ν), µ = νθ, (18)
with

ν = 1 + 3S + 9S 2 + 16(S + 1) 4 , S = G K . (19)
Moreover, we have the following relation: µ

+ θ 2 = K + G.
Remark 1. We highlight the fact that F → ⟨Q BGK ⟩ [F ] = ⟨M⟩ [F ] -F and the relaxation operator F → G[F ] -F are two different operators. As mentioned in [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF], the average operator ⟨Q BGK ⟩ is more complicated than a relaxation type operator defined with gyromoments.

Remark 2. Let us remark that P := x + ⊥ v is invariant under the change of variable i.e. x

+ ⊥ v = X + ⊥ V .
Thus, in this set of variables (P , V ), we have dx dv = dX dV = dP dV and the above calculations can be performed by replacing F (X, V ) by F (P , X 3 , V ) in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF] and the gyromaxwellian can be written in a more natural way by

G[F ](P , X 3 , V ) = nM 2 µθ µ-θ (V -U ) M 1 θ (V 3 -U 3 ) M 2 µ (P -Y ). ( 20
)
The variables (P , X 3 , V ) turn out to be natural for the gyrokinetic model [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF].

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

Discretisation of models

In this section, we develop multiscale numerical schemes for the kinetic model [START_REF] Birdsall | Plasma physics via computer simulation[END_REF] or equivalently [START_REF] Bostan | High magnetic field equilibria for the Fokker-Planck-Landau equation[END_REF]. The goal is to design numerical schemes which overcome the stiffness induced by the two physical parameters ε and τ . 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 ε using scales-separation strategy developed in [5-9, 14, 15] which introduces a new time periodic variable; the dissipative stiffness in τ is handled using IMEX and exponential time integrations which are widely used to achieve asymptotic preserving properties [START_REF] Dimarco | Exponential Runge-Kutta Methods for Stiff Kinetic Equations[END_REF][START_REF] Hochbruck | Exponential integrators[END_REF][START_REF] Jin | A Successive Penalty-Based Asymptotic-Preserving Scheme for Kinetic Equations[END_REF].

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 τ and ε to alleviate notations. Let N p ∈ N a number of numerical macro-particles of positions x q , velocities v q and weights ω q , q = 1, . . . , N p . We will approximate the solution f of (2) by the following particle ansatz

f Np (t, x, v) = Np q=1 ω q (t)δ(x -x q (t))δ(v -v q (t)) = Np q=1 ω q (t)δ(z -z q (t)), (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 caracteristics of [START_REF] Bostan | High magnetic field equilibria for the Fokker-Planck-Landau equation[END_REF].

żp (t) = 1 ε Az p (t) + h(t, z p (t)), z p (0) = z p,0 . (22)
We assume that z p,0 ∈ Ω ∪ D v ⊂ R 6 where Ω ⊂ R 3 is the spatial domain and D v ⊂ R 3 the numerical velocity domain. Due to the presence of the source term

Q[f ](t, z) = 1 τ (M[f ](t, z) -f (t, z))
, the weights ω q (t) in [START_REF] Gamba | Spectral -Lagrangian methods for Collisional Models of Non -Equilibrium Statistical States[END_REF] are not constant in time. We denote by s p (t) = 1 τ (m p (t) -ω p (t)) the weights of the source term with the Maxwellian's weights

m p (t) = M[f ](t, z p (t))V Np and ω p (t) = f (t, z p (t))V Np the weights of f , with V Np = |Ω||D v |/N p
is the average volume of the phase space domain occupied by one particle. To compute the moments of f Np , we need to regularize f Np with respect to x. In this work we consider tensorial product of first order B-splines: 3 and N x the grid point number, we choose ∆x = (x max -x min )/N x , then we have

S(x 1 , x 2 , x 3 ) = S(x 1 ) S(x 2 ) S(x 3 ). If Ω = [x min , x max ]
S(x) =      1 ∆x 1 - |x| ∆x if |x| ≤ ∆x, 0 else. ( 23
)
With the smooth kernel S, we define the regularized form

f S Np of f Np f S Np (t, x, v) = Np q=1 ω q (t)S(x -x q (t))δ(v -v q (t)), (24)
from which we can compute the moments of f S Np following n, nu, n |u|

2 2 + 3 2 T (t, x) = Np q=1 ω q (t)S(x -x q (t))(1, v q (t), |v q (t)| 2 /2), (25)
and then, we can define M[f S Np ] (using ( 11)) together with its weights

m p (t) = M[f S Np ](t, z p (t))V Np . (26)
Notice that m p (t) does not only depend on z p , but also depends on all weights and particles through f S Np 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 ω p , which determine f S Np ). Finally, the weights ω p satisfy the following equation

ωp (t) = 1 τ (m p (t) -ω p (t)). ( 27
)
Now if we replace collisional operator by its particles version 1 τ q (m q -ω q )δ(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 [START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF][START_REF] Issautier | Convergence of a weighted particle method for solving the boltzmann (bgk) equation[END_REF][START_REF] Sonnendrücker | Numerical Methods for the Vlasov-Maxwell equations[END_REF]). The goal is then to solve numerically the multiscale system ( 22)-( 27) endowed with initial conditions (z p (0), ω p (0)) = (z p,0 , ω p,0 ).

As in the continuous case, the main oscillation generated by żp = (1/ε)Az p is filtered out by a change of unknown. We set Z p (t) = e -t ε A z p (t) (with A given by ( 3)) so that [START_REF] Hochbruck | Exponential integrators[END_REF] becomes

Żp (t) = H(t, t/ε, Z p (t)) := e -t/εA h(t, e t/εA Z p (t)). (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 Np (t, Z = (X, V )) = Np q=1 ω q (t)δ(X -X q (t))δ(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 S Np of F Np thanks to the unfiltered particle density f S Np with filtered variables:

F S Np (t, Z) = f S Np (t, X (t/ε, Z) , V (t/ε, Z)) = Np q=1 ω q (t)S (X (t/ε, Z)-X (t/ε, Z q (t))) δ (V (t/ε, Z)-V (t/ε, Z q (t))) . (29)
We are now able to compute moments of F S Np as

N (t, Z) = Np q=1 ω q (t)S(X(t/ε, Z) -X(t/ε, Z q (t))), N U(t, Z) = Np q=1 ω q (t)S(X(t/ε, Z) -X(t/ε, Z q (t)))V(t/ε, Z q (t)), N |U| 2 2 + 3 2 T (t, Z) = Np q=1 ω q (t)S(X(t/ε, Z) -X(t/ε, Z q (t))) |V(t/ε, Z q (t))| 2 2 , (30)
and to reconstruct the corresponding filtered Maxwellian M filt [F S Np ]. Let us remark that we have (N , U, T )(t, Z) = (n, u, T )(t, X(t/ε, Z)) and that the weights of the Maxwellian M and the filtered Maxwellian M filt [F ] are equal since the filtration only applies on the particles

M[f S Np ](t, z) = Np q=1 m p (t)δ(z -z q (t)), M filt [F S Np ](t, Z) = Np q=1 m p (t)δ(Z -Z q (t)),
and the same holds true for f Np and F Np . Finally, once the filtered Maxwellian M filt [F S Np ](t, Z) has been constructed, the weights are computed by evaluating it at the particles Z p (t):

m p (t) = M filt [F S Np ](t, Z p (t))V Np .
To anticipate with the next subsection devoted to the numerical schemes, we introduce a new definition to highlight the "operational" dependence of M filt (with respect to the weights and particles positions) and its dependence with respect to the fast variable

m p (t) = M filt [Z(t), ω(t)](t, t/ε, Z p (t))V Np , (31)
where Z(t) = {Z q (t)} q=1,...,Np and ω(t) = {ω q (t)} q=1,...,Np .

Multiscale numerical schemes

The presence of a slow time scale and a fast one which is 2π-periodic is amenable to the two-scale framework [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Chartier | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF][START_REF] Chartier | Uniformly accurate methods for three dimensional Vlasov equations under strong magnetic field with varying direction[END_REF][START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF][START_REF] Chartier | Derivative-free high-order uniformly accurate schemes for highly-oscillatory systems[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF]. Then, we consider the two time scales as independent by introducing Z p (t, s)

(resp. W p (t, s)) satisfying Z p (t, t/ε) = Z p (t) (resp. W p (t, t/ε) = ω p (t)
). We also introduce the following functions:

H p (t, s) = H(t, s, Z p (t, s)), M p (t, s) = M filt [Z(t, s), W(t, s)](t, s, Z p (t, s))V Np . (32)
With these notations, ( 28) and ( 27) become

∂ t Z p (t, s) + 1 ε ∂ s Z p (t, s) = H p (t, s), ∂ t W p (t, s) + 1 ε ∂ s W p (t, s) = 1 τ (M p (t, s) -W p (t, s)), (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 [START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF]. In this part, we present different multiscale numerical schemes according to the considered asymptotic in ε or/and τ . 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π/N s , j = 0, . . . , N s , with N s ∈ N ⋆ (an even number), the Fourier transform writes, for an arbitrary function g : s ∈ [0, 2π] → g(s) of Fourier coefficients ĝj ,

g(s) = Ns/2-1 l=-Ns/2 ĝl e ils .
Applying the Fourier transform to [START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF] enables to get the following equations on the modes

d dt Ẑp,l (t) + il ε Ẑp,l (t) = Ĥp,l (t), d dt Ŵp,l (t) + il ε + 1 τ Ŵp,l (t) = 1 τ Mp,l (t).
We introduce time discretisation t n = n∆t, ∀n ∈ N with ∆t > 0 the time step. Applying the variation of constant formula between t n and t n+1 leads to

Ẑp,l (t n+1 ) = e -il ε ∆t Ẑp,l (t n ) + t n+1 t n e -il ε (t n+1 -t) Ĥp,l (t) dt, Ŵp,l (t n+1 ) = e -( il ε + 1 τ )∆t Ŵp,l (t n ) + 1 τ t n+1 t n e -( il ε + 1 τ )(t n+1 -t)
Mp,l (t) dt.

(34)

These expressions will pave the way of the time integrators that will be described in the sequel.

Scheme A: uniform accuracy with respect to ε

In this subsection, we introduce a first numerical scheme to approximate [START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF]. The use of the two-scale approach will make this scheme uniformly accurate with respect to ε for all fixed τ > 0. Besides, this scheme enjoys the asymptotic preserving property in the asymptotic τ → 0 for all fixed ε > 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 [START_REF] Chartier | Derivative-free high-order uniformly accurate schemes for highly-oscillatory systems[END_REF], [START_REF] Mocquard | jl: A Julia package for highly oscillatory problems[END_REF]. Denoting by g(t) a function representing Ĥp,l (t) or Mp,l (t) in (34), we replace g(t) by the following approximations g(t) ≈ g(t n ) (first order approximation),

g(t) ≈ g(t n ) + (t -t n ) g(t n ) -g(t n-1
) ∆t (second order approximation).

To simplify the final expression of the scheme, we introduce for all α ∈ N, z ∈ C, the following function

φ α (z) := 1 α! 1 0 e -z(1-t) t α dt = (-1) α+1 z α+1 e -z - α k=0 (-z) k k! , so that t n+1 t n e -z(t n+1 -t) (t -t n ) α dt = α! ∆t α+1 φ α (z∆t). ( 35
)
We

Ŵn+1 p,l = e -( il ε + 1 τ )∆t Ŵn p,l + ∆t τ φ 0 il ε + 1 τ ∆t Mn p,l , (36a)
and the second order one is

Ẑn+1 p,l = e -il ε ∆t Ẑn p,l + ∆tφ 0 ∆t il ε Ĥn p,l + ∆tφ 1 ∆t il ε ( Ĥn p,l -Ĥn-1 p,l ), Ŵn+1 p,l = e -( il ε + 1 τ )∆t Ŵn p,l + ∆t τ φ 0 ∆t il ε + 1 τ Mn p,l + ∆t τ φ 1 ∆t il ε + 1 τ ( Mn p,l -Mn-1 p,l ). ( 37
)
Then, once the two-scale (Fourier) numerical solutions are computed, one can evaluate the Fourier series

Z n p (s) = Ns/2-1 l=-Ns/2 Ẑn p,l e ils , W n p (s) = Ns/2-1 l=-Ns/2 Ŵn p,l e ils , (38)
and the original solutions are obtained by evaluating on the diagonal s = t n /ε and applying the inverse change of variable

Z n p = Z n p t n ε , (x n p , v n p ) = z n p = e t n ε A Z n p , ω n p = W n p t n ε . ( 39
)
Let us discuss the initial condition. By introducing the variable s, we obtain a degree of freedom on the initial condition for [START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF]. The only constraint is Z p (0, 0) = Z p (0) and W p (0, 0) = ω p (0). We use this degree of freedom to control the smoothness of Z p (t) and W p (t). To do so, well-prepared initial conditions can be constructed that allow the scheme to be uniformly accurate with respect to ε. For the first order approximation, we consider constant in s initial condition

Z p (0, s) = Z p (0) = (x p (0), v p (0)), W p (0, s) = ω p (0). ( 40
)
For the second order approximation, we consider

Z p (0, s) = Z p,0 + ε s 0 (I -Π)H(0, θ, Z p,0 ) dθ, W p (0, s) = ω p,0 + ε τ s 0 (I -Π)M p (0, θ) dθ, (41)
where Πg = 1 2π 2π 0 g(s)ds is the averaging operator (see [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF] for more details). This preparation is only required to ensure uniform estimates of the time derivatives with respect to ε and is not used in the classical fluid limit (ε > 0 fixed and τ → 0). The two-scale framework provides uniform accuracy for (36) and (37) with respect to ε, for fixed τ . Conversely, for a fixed ε and looking for the limit τ → 0, the problem becomes nonstiff in ε and dissipative in τ . For such a problem, several numerical schemes have been proposed, entering in the category of Asymptotic Preserving schemes [START_REF] Jin | Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations[END_REF]. This is the purpose of the next part.

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)

Ŵn+1 p,l = e -( il ε + 1 τ )∆t Ŵn p,l + ε ilτ + ε 1 -e -( il ε + 1 τ )∆t M⋆ p,l ,
with M⋆ p,l = Mn p,l . We have an AP property in the following sense: when τ → 0, we have Ŵn+1

p,l
→ M⋆ p,l and thus our density weights converge toward maxwellian weights. However to get a better asymptotic behavior when τ → 0, an implicit scheme is required by taking M⋆ p,l = Mn+1 p,l . However, computing Mn+1 p,l requires the knowledge of Ẑn+1 p,l and Ŵn+1 p,l according to [START_REF] Sonnendrücker | Numerical Methods for the Vlasov-Maxwell equations[END_REF]. Actually, Mp,l depends on the weights through the moments, but here, contrary to the usual Maxwellian case [START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Jin | Efficient Asymptotic-Preserving (AP) Schemes For Some Multiscale Kinetic Equations[END_REF], the moments [START_REF] Ren | An asymptotic-preserving Monte Carlo method for the Boltzmann equation[END_REF] 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 ε. Thus, to avoid a costly fully implicit scheme, we consider a hybrid approximation of (31) computed from the knowledge of Ẑn+1 p,l and Ŵn p,l . The procedure is then the following: we first compute Z n+1 q (s) from the first equation of (36) and then we compute the moments according to

U ⋆ (s, x) = Np q=1 W n p (s)S(x -X(s + ∆t ε , Z n+1 q (s + ∆t ε )))m(V(s + ∆t ε , Z n+1 q (s + ∆t ε ))), (42)
where m(v) = (1, v, |v| 2 /2). Then, we are able to compute the Maxwellian weights M ⋆ p (s) according to these moments:

M ⋆ p (s) = M filt [Z n+1 (s + ∆t/ε), W n (s)](t n+1 , s + ∆t/ε, Z n+1 p (s + ∆t/ε))V Np . ( 43
)
Once the weights of the filtered Maxwellien are computed, we then update the weights with

Ŵn+1 p,l = e -( il ε + 1 τ )∆t Ŵn p,l + ε ilτ + ε 1 -e -( il ε + 1 τ )∆t M⋆ p,l , (44)
with M⋆ p,l computed with (43). Then, when τ → 0, Ŵn+1 p,l converges towards M⋆ p,l . For the second order scheme, when τ goes to 0, Ŵn+1 p,l converges toward 2 Mn p,l -Mn-1 p,l , which is the linear extrapolation of Mn+1 p,l .

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 ⟨Q BGK ⟩. This important property is not ensured by the scheme A presented in Subsections 3.3 and 3.4. When ε → 0, the averaged BGK operator ⟨Q BGK ⟩ [F ] = ⟨M filt ⟩[F ] -F turns out to be more complicated than standard BGK-like operator. Indeed, as τ → 0, F converges towards G[F ] and not

⟨M filt ⟩[F ].
Hence, the averaged BGK operator ⟨Q BGK ⟩ shares similarities with more complex nonlinear operators like the Boltzmann or Landau operator for which specific Asymptotic Preserving schemes have been developed (see [START_REF] Dimarco | Exponential Runge-Kutta Methods for Stiff Kinetic Equations[END_REF][START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]) 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 G[F ] defined in [START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF]. In this section, we focus on the relaxation part (τ → 0) and thus, we only consider the equation on the weights, the particles positions being approximated by (36a). Starting from the equation [START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF] on the weights, we add and substract the gyromaxwellian weights (denoted by (1/τ ) G p ). We compute G p (t, s) in a similar way than maxwellian weights: we compute gyromoments associated to W p , which determines a unique gyromaxwellian G. Then we compute gyromaxwellian weights as in [START_REF] Wang | Geometric two-scale integrators for highly oscillatory system: uniform accuracy and near conservations[END_REF] with G instead of M filt . Equation [START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF] on the weights becomes

∂ t W p + 1 ε ∂ s W p = 1 τ (M p -W p ) - 1 τ G p + 1 τ G p . ( 45
)
Considering the Fourier transform in s, the variation of constant formula leads to

∂ t Ŵp,l e ( il ε + 1 τ )t = e ( il ε + 1 τ )t 1 τ ( Mp,l -Ĝp,l ) + 1 τ Ĝp,l . (46)
We now integrate between t n and t n+1 to get

Ŵp,l (t n+1 ) = e -( il ε + 1 τ )∆t Ŵp,l (t n ) + 1 τ t n+1 t n e -( il ε + 1 τ )(t n+1 -t) ( Mp,l (t) -Ĝp,l (t)) dt + 1 τ t n+1 t n e -( il ε + 1 τ )(t n+1 -t) Ĝp,l (
To compute Ĝn+1 p,l , we need to compute the gyromoments U n+1 defined by [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF]. 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

∂ t U + 1 ε ∂ s U = 1 τ ( Ũ -U),
where Ũ denotes the gyromoments of the filtered Maxwellian M filt ( whose weights are denoted by M p ). An exponential scheme to update the Fourier modes of U from t n to t n+1 is given by

Ûn+1 l = e -( il ε + 1 τ )∆t Ûn l + ε ilτ + ε (1 -e -( il ε + 1 τ )∆t ) Ûn l , l ̸ = 0, Ûn+1 0 = e -∆t τ Ûn 0 + (1 -e -∆t τ ) Ûn 0 . (48)
As we will see in the next section, this scheme is AP in the gyrokinetic limit (ε → 0 and τ > 0 fixed) and in the gyrofluid limit (ε → 0 then τ → 0). However, for this scheme, the classical fluid limit (τ → 0 and ε > 0 fixed) is not recovered since ∀ε ≥ 0, Ŵn+1 p,l → Ĝn+1 p,l when τ → 0.

Convex combined scheme

The scheme B presented in the last subsection does not capture the correct limit when τ → 0 for a fixed ε since it converges to a gyromaxwellian by construction. With the notation Ŵn+1,A p,l obtained from (44) and Ŵn+1,B p,l obtained from (47)-( 48), we consider a convex combination of these two numerical schemes with coefficients ε τ +ε and 1 -ε τ +ε so that, according to the desired regimes, the coefficients degenerate towards the correct scheme:

Ŵn+1 p,l = ε τ + ε Ŵn+1,A p,l + 1 - ε τ + ε Ŵn+1,B p,l .
To summarise the different asymptotics of this convex combined scheme, when ε > 0 fixed and τ → 0, the scheme degenerates into scheme A which has right behaviour in the classical fluid limit. When τ > ε → 0, the scheme degenerates into scheme B which has the good behaviour in the gyrofluid limit.

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 (ε → 0 then τ → 0). In the second subsection, we show that, when the Maxwellian is independent of weights, scheme A is UA in ε and in τ .

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 ε → 0 then τ → 0.

Proof. We first consider the limit ε → 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 ε → 0 except for the mode l = 0). Then, in (36a), we obtain when ε → 0 Ẑn+1 p,0 = Ẑn p,0 + ∆t Ĥn p,0 , whereas for the weights equations ( 47

Ûn+1 0 = e -∆t τ Ûn 0 + (1 -e -∆t τ ) Ûn 0 = Ûn 0 , (49)
since Ûn 0 = Ûn 0 when ε = 0 (indeed, ⟨Q BGK ⟩ conserves gyromoments). This scheme is a consistent explicit first order time discretization of the gyrokinetic equation [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF].

We now focus on the limit τ → 0. Following [START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Ren | An asymptotic-preserving Monte Carlo method for the Boltzmann equation[END_REF][START_REF] Zhang | An asymptotic preserving Monte Carlo method for the multispecies Boltzmann equation[END_REF] 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 τ → 0 towards the gyromaxwellian weights. It is clear from (49)-(50) that Ŵn+1 p,0 converges to Ĝn+1 p,0 which means that the asymptotic scheme is a consistent discretization of the gyrofluid equations (65) (presented in Appendix A).

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 Z only. In this case, one can prove uniform accuracy in time with respect to both ε and τ . When the Maxwellian depends on the weights W, we are able to prove uniform accuracy in time with respect to ε (using classical tools from [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF], [START_REF] Wang | Geometric two-scale integrators for highly oscillatory system: uniform accuracy and near conservations[END_REF]) and an asymptotic property in τ .

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

∥∂ k t Z p (t, •)∥ L ∞ s ≤ C, k = 0, 1, 2. ( 51
)
We consider the model for (Z p , W p ) : (t, s) ∈ [0,

t max ] × [0, 2π] → (Z, W)(t, s) ∈ R 6Np × R Np , ∂ t Z p (t, s) + 1 ε ∂ s Z p (t, s) = H p (t, s), ∂ t W p (t, s) + 1 ε ∂ s W p (t, s) = 1 τ (M p (t, s) -W p (t, s)), (52)
where H p (t, s) = H(t, s, Z p (t, s)) ∈ R 6Np and M p (t, s) = M (t, s, Z p (t, s)) ∈ R Np , with a well-prepared initial condition Z p (0, s), W p (0, s) computed as in (41). Obviously, some assumptions are required. In the sequel, the norm ∥ • ∥ of a vector refers to the standard euclidean norm in R d (d ∈ 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 ω. 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 ε and τ : there exists a constant C > 0, independent of ε and τ such that

∥Z n (t n /ε) -z(t n )∥ L ∞+ ∥W n (t n /ε) -w(t n )∥ L ∞ ≤ C∆t 2 , 0 ≤ n ≤ N, N ∆t = t max , ( 53
)
where z(t), w(t) are the solutions of ( 22) and [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF].

Proof. Recall that H p (t, s) = H(t, s, Z p (t, s)), M p (t, s) = M filt [Z(t, s), W(t, s)](t, s, Z p (t, s)), and
H n p (s) = H(t n , s, Z n p (s)), M n p (s) = M filt [Z n (s), W n (s)](t n , s, Z n p (s)
). The solution Z p of the two-scale equation (52) satisfies

Z p (t n+1 , s) = l∈Z Ẑp,l (t n+1 )e ils = l∈Z   e -il ε ∆t Ẑp,l (t n ) + t n+1 t n e -il ε (t n+1 -t) Ĥp,l (t) dt    e ils . (54)
A second order Taylor-Lagrange formula enables to expand the term Ĥp,l (t) at t = t n to get

Z p (t n+1 , s) = l∈Z   e -il ε ∆t Ẑp,l (t n ) + t n+1 t n e -il ε (t n+1 -t) Ĥp,l (t n ) dt + t n+1 t n e -il ε (t n+1 -t) (t -t n ) d dt Ĥp,l (t n ) dt    e ils + l∈Z    t n+1 t n e -il ε (t n+1 -t) (t -t n ) 2 1 0 (1 -θ) d 2 dt 2 Ĥp,l (t n + θ(t -t n )) dθ dt    e ils .
Using l∈Z e ils Ĥp,l (t n ) = H p (t n , s) and a Taylor-Lagrange formula on d dt Ĥp,l (t n ), we get

Z p (t n+1 , s) = Z p (t n , s - ∆t ε ) + t n+1 t n H p (t n , s - t n+1 -t ε ) dt + t n+1 t n (t -t n ) ∆t H p (t n , s - (t n+1 -t) ε ) -H p (t n-1 , s - (t n+1 -t) ε ) dt + r n p (s), ( 55
)
where r n p (s) denotes the truncature error defined by

r n p (s) = l∈Z    t n+1 t n e -il ε (t n+1 -t) (t -t n ) 1 0 (1 -θ)∆t d 2 dt 2 Ĥp,l (t n -θ∆t) dθ dt    e ils + l∈Z    t n+1 t n e -il ε (t n+1 -t) (t -t n ) 2 1 0 (1 -θ) d 2 dt 2 Ĥp,l (t n + θ(t -t n )) dθ dt    e ils = t n+1 t n 1 0 (1 -θ)(t -t n )∆t d 2 dt 2 H p (t n -θ∆t, s - t n+1 -t ε ) dθ dt + t n+1 t n 1 0 (1 -θ)(t -t n ) 2 d 2 dt 2 H p (t n + θ(t -t n ), s - t n+1 -t ε ) dθ dt. ( 56
)
Let us recall that the time derivative of Ĥp,l also involves time derivatives of Z p and derivatives ∂ z Ĥp,l . With the assumptions performed above, these quantities are all uniformly bounded. In the same way, the numerical scheme (37) on Z rewrites

Z n+1 p (s) = Z n p s - ∆t ε + l∈Z ∆tφ 0 il ε ∆t Ĥn p,l +∆t 2 φ 1 il ε ∆t Ĥn p,l -Ĥn-1 p,l ∆t = Z n p s - ∆t ε + t n+1 t n H n p s - t n+1 -t ε dt + t n+1 t n (t -t n ) ∆t H n p s - t n+1 -t ε -H n-1 s - t n+1 -t ε dt, (57)
where we recall that ∆t α+1 φ α il ε ∆t =

t n+1 t n e -il ε (t n+1 -t) (t -t n ) α dt for α = 0, 1.
Now, in view of deriving an induction on the error e n Zp (s) = Z p (t n , s) -Z n p (s), we define, for n ≤ N

∆H n p (s) = H p t n , s - t n+1 -t ε -H n p s - t n+1 -t ε .
With this, considering the difference between (55) and (57), we obtain the following induction on e n Zp (s)

e n+1 Zp (s) = e n Zp s - ∆t ε + t n+1 t n ∆H n p (s) dt + t n+1 t n (t -t n ) ∆t ∆H n p (s) -∆H n-1 p (s - ∆t ε ) dt + r n p (s). ( 58
)
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 ε and τ because of (51). Under the assumptions on the nonlinear term H, we have

∆H n p (s) ≤ C e n Zp (s) ≤ C e n Zp L ∞ .
Moreover, the uniform estimates of ∂ k t Z p (k = 0, 1, 2) ensure that the local error can be estimated by

r n p (s) ≤ ∆t 3 d 2 dt 2 H p L ∞ ≤ C∆t 3 . ( 59
)
Setting t max = N ∆t and taking the maximum on p leads to

e n+1 Z L ∞ ≤ (1 + C∆t)∥e n Z ∥ L ∞ + C∆t e n-1 Z L ∞ + C∆t 3 . ( 60
) Let us denote U n = (∥e n Z ∥ L ∞ , e n-1 Z L ∞ ) T , A = 1 + C∆t C∆t 1 0 , b = C∆t 3 0 .
We write (60) as

U n+1 ≤ AU n + b to get recursively U n ≤ A n U 0 + n-1 k=0 A k b.
Taking the norm, we get

U N ≤ (1 + 2C∆t) N U 0 + N -1 k=0 (1 + 2C∆t) k b.
Once again, we denote by C the constant involved. Since U 0 = max(∥e 0 ∥, ∥e -1 ∥) = ∥e -1 ∥ = C∆t 2 (iteration n = -1 computed with the order 1 scheme), we have finally

U N ≤ Ce CT ∆t 2 + e CT N -1 k=0 b ≤ Ce CT ∆t 2 .
Let consider now the numerical scheme on the weights. First, using the same calculations (Taylor-Lagrange formula on Mp,l (t)) as before, the two-scale solution W p of (52) satisfies

W p (t n+1 , s) = e -∆t τ W p (t n , s - ∆t ε ) + 1 τ l∈Z    t n+1 t n e -( il ε + 1 τ )(t n+1 -t) Mp,l (t n ) dt    e ils + 1 τ l∈Z    t n+1 t n e -( il ε + 1 τ )(t n+1 -t) (t -t n ) Mp,l (t n ) -Mp,l (t n-1 ) ∆t dt    + q n p (s), (61)
where q n p (s) denotes the truncature error

q n p (s) = 1 τ l∈Z    t n+1 t n e -( il ε + 1 τ )(t n+1 -t) (t -t n ) 1 0 (1 -θ)∆t d 2 dt 2 Mp,l (t n -θ∆t) dθ dt   e ils + 1 τ l∈Z    t n+1 t n e -( il ε + 1 τ )(t n+1 -t) (t -t n ) 2 1 0 (1 -θ) d 2 dt 2 Mp,l (t n + θ(t -t n )) dθ dt    e ils . (62)
Considering the difference between (61) and the scheme (37) on W leads to an induction on the error

e n+1 Wp (s) = e -∆t τ e n Wp s - ∆t ε + 1 τ t n+1 t n e -t n+1 -t τ ∆M n p (s) dt + 1 τ t n+1 t n e -t n+1 -t τ t -t n ∆t dt ∆M n p (s) -∆M n-1 p (s) + q n p (s), (63)
where we introduced the following notation for n ≤ N

∆M n p (s) = M p t n , s - t n+1 -t ε -M n p s - t n+1 -t ε .
Thanks to the definition of the φ function (35) and the assumptions on M , the truncature error q n p (s) given by (62) satisfies the following estimate

q n p (s) ≤ d 2 M p dt 2 L ∞ 1 τ    t n+1 t n e -1 τ (t n+1 -t) (t -t n ) ∆t dt + t n+1 t n e -1 τ (t n+1 -t) (t -t n ) 2 dt    ≤ d 2 M p dt 2 L ∞ ∆t 3 τ φ 1 ∆t τ + 2φ 2 ∆t τ .
Consequently, under the assumptions on M and the uniform estimates on e Zp , we have ∆M n p (s) ≤ C∥e n Z ∥ L ∞ ≤ C∆t 2 . Then, the second term in the right hand side of (63) can be estimated by

1 τ t n+1 t n e -t n+1 -t τ ∆M n p (s) dt ≤ C ∆t 3 τ φ 0 ∆t τ ,
whereas the third term satisfies

1 τ t n+1 t n e -t n+1 -t τ (t -t n ) dt ∆M n p (s) -∆M n-1 p (s) ∆t ≤ 2C ∆t 3 τ φ 1 ∆t τ .
Thus, (63) becomes

e n+1 W L ∞ ≤ e -∆t τ ∥e n W ∥ L ∞ + C ∆t 3 τ φ 0 ∆t τ + 3φ 1 ∆t τ + 2φ 2 ∆t τ .
By induction, we have, with

N ∆t = t max , e N W L ∞ ≤ e -tmax/τ e 0 W L ∞ + C ∆t 3 1 -e -tmax/τ τ 1 -e -∆t/τ (φ 0 + 3φ 1 + 2φ 2 ) ∆t τ , ≤ e -tmax/τ e 0 W L ∞ + C(1 -e -tmax/τ ) ∆t 2 +2τ 2 -3τ ∆t+ 4∆t 2 -2τ ∆t 1 -e -∆t/τ .
One can prove that z 2 + 2 -3z + 4z 2 -2z 1-e -z ≤ 5z 2 for all z > 0 from which we deduce, using e 0 W L ∞ = 0,

e N W L ∞ ≤ e -tmax/τ e 0 W L ∞ + C(1 -e -tmax/τ )5∆t 2 ≤ C(1 -e -tmax/τ )∆t 2 ,
which concludes the proof.

Remark 3. When M depends on ω, in the truncature error (62) the term d 2 dt 2 M makes appear ∂ 2 t W. For a fixed τ , time derivatives are uniformly bounded in ε, which enable to garantee the UA property in ε. However, when τ → 0 the time derivatives of W are not uniformly bounded in τ . In this case, we thus lose uniform accuracy in τ on weights (this property is still present on particles since it does not depend on τ ). For the AP property when τ → 0, we refer to Subsection 3.4.

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.

ODE framework

We consider the following multiscale ODE satisfied by z = (x, v) ∈ R 6 and ω ∈ R

ż(t) = 1 ε Az(t) + h(t, z), ω(t) = 1 τ (M(ω, z) -ω(t)), (64)
with A and h given by ( 3) and E(t, x) = ((x 1 + x 3 ) cos(t), x 1 x 2 sin(t), -x 2 2 e -t 2 ). Two examples for M will be investigated. In a first example, M(ω, z) = 1 + |z| 2 + e -ω 2 (M depends on ω) whereas in a second example

M(z) = 1 + |z| 2 (M does not depend on ω). The initial condition is z(0) = (1, 1, 0, 1/2, 1/2, 3/2)
, ω(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 ∆t = 2 -14 . The error considered is the sum of the L ∞ -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 M(ω, z). We plot the error as a function of the time step ∆t = 2 -k , k = 3, . . . , 13 for different values of ε = 10 -k (k = 0, . . . , 6) with τ = 1 (left figure) and of τ = 10 -k (k = 0, . . . , 4) with ε = 1 (right figure). We observe the uniform accuracy with respect to ε (left figure) since for all the values of ε, the error remains the same and the second order is recovered, as expected. However, it is not the case when τ varies (right figure) since the order is observed only when τ > ∆t. In Figure 2, we consider the case M(z) for which Theorem 1 ensures uniform accuracy in both ε and τ parameters. As expected, we obtain the uniform accuracy property in τ (the uniform accuracy in ε is ensured by the two-scale strategy). L -error Z + W eps=10 0.0 eps=10 1.0 eps=10 2.0 eps=10 3.0 eps=10 4.0 eps=10 5.0 eps=10 6.0 (a) L ∞ -Error for τ = 1 and different ε. In Figures 3, we plot the time history of the Maxwellian weights M(z(t)) together with w(t) for different values of ε and τ , with ∆t = 2.5 × 10 -4 . These two figures illustrate the asymptotic properties (gyrokinetic limit ε → 0 and fluid limit τ → 0) of (64). Indeed, when ε is fixed, we observe the convergence when τ → 0 of w(t) to M(z(t)) which oscillates in time (see Figure 3-(a)). However, for τ > 0 fixed, and ε → 0, we observe the convergence of w(t) towards the average of M(z(t)). Finally, in Figure 4, we illustrate the Uniform Accuracy property in both parameters ε and τ . We take ε = τ = 0.01 and plot the time history of ω for two different values of the time step: ∆t = 2.5 × 10 -4 and ∆t = 0.05. We can observe that even if ∆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.

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 ε,τ (t, Z) becomes

∂ t F ε,τ = 1 τ Q filt [F ε,τ ](t, t/ε, Z) = 1 τ (M filt [F ε,τ ] -F ε,τ ),
where M filt is defined in [START_REF] Crestetto | A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling[END_REF]. 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 (ε → 0 then τ → 0). Hence, from an out-of-equilibria initial data F 0 , we compute its associated equilibria G[F 0 ], and observe the relaxation of F ε,τ towards 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 , Y a , K a , G a ) and

U b = (n, U b , Y b , K b , G b ).
After some computations, we find that the gyromaxwellian G[F 0 ] associated to F 0 is characterized by

U ab = 2n, Ua+U b 2 , Y a+Y b 2 , Ka+K b 2 + |Ua-U b | 2 8 , Ga+G b 2 + |Y a -Y b | 2 8 - |Ua-Ub| 2 8
.

For the first test, we consider the scheme B (36a)-( 47)-(48) with the following numerical parameters: N p = 1.6×10 7 particles, N s = 8 for the s direction and ∆t = 2 × 10 -4 whereas ε = 10 -8 and τ = 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 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 ε,τ (t, Z), we plot a reconstruction of F ε,τ (t = 0, X 1 , X 2 = 0, V 1 = 0, V 2 = 0, V 3 = 0) and F ε,τ (t = t 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 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. 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 n 1 = p |ω n p -G ex p |, e n 2 = p |ω n p -G n p | and e n 3 = p |G n p -G ex p |
, where G ex p denotes the weights computed from analytical gyromaxwellian G[F 0 ] (which are constant in time) and G n p 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, ε = 10 -8 , ∆t = 10 -3 , N s = 4, ∆x ≈ 0.37 and we look at different values of τ .

In Figures 6 and7, we can see that e n 3 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 G n p . 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 [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF][START_REF] Gamba | Spectral -Lagrangian methods for Collisional Models of Non -Equilibrium Statistical States[END_REF] to enforce the conservation of the gyromoments at the discrete level.

Let us now focus on the other errors e n 1 and e n 2 for the two schemes. In Figure 6, the error e n 2 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 ∆t > τ ) thanks to its AP property. The error e n 1 is the error e n 2 with the effect of non conservation of gyromoments which are deteriorated with time.

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 [START_REF] Bostan | The finite Larmor radius regime : collisional setting and fluid models[END_REF], involves two scaling parameters: ε drives the oscillation due to the presence of the strong external magnetic field and τ 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 [START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. Considering the coupling with the self-consistent electromagnetic field governed by Maxwell equation could be also considered.

A Gyrofluid model's properties

In this section, we study the gyrofluid model derived from the gyrokinetic model ( 8) with the BGK operator when τ → 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 τ by the gyromaxwellian (15) enables to get the following gyrofluid model satisfied by U(t,

X 3 ) = (n, nU , nU 3 , nY , n(K + |U | 2 /2), n(G + |Y | 2 /2 -|U | 2 /2))(t, X 3 ) ∈ R 8 ∂ t n + ∂ X3 nU 3 = 0, (65a)
∂ t nU + ∂ X3 nU 3 U = 0, (65b)
∂ t nU 3 + ∂ X3 nU 2 3 + nθ = nE 3 , (65c)
∂ t nY + ∂ X3 nU 3 Y = n ⊥ E, (65d)
∂ t n K + |U | 2 2 + ∂ X3 nU 3 K + |U | 2 2 + θ = nU 3 E 3 , (65e)
∂ t n G + Y 2 2 - U 2 2 + ∂ X3 nU 3 G + Y 2 2 - U 2 2 = nY • ⊥ E. (65f)
The gyrofluid model has also an entropic inequality, derived directly from H-theorem's inequality

∂ t n ln n µ 2 θ 3/2 µ-θ + ∂ X3 nU 3 ln n µ 2 θ 3/2 µ-θ ≤ 0, (66)
µ and θ being defined by [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. 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. ( Y 2 -U 2 )/2). Moreover, as θ is a function of the two internal energies K and G, nθ 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, Y , K, G)(t, X 3 ) ∈ R 8 . This is done in the following Proposition.

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

∂ t W + A(W)∂ X3 W = B, (67)
with W(t, X 3 ) = (n, U, Y , K, G)(t, X 3 ) and with

A(W) =             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 θ n 0 0 U 3 0 0 θ K θ G 0 0 0 0 U 3 0 0 0 0 0 0 0 0 U 3 0 0 0 0 0 θ 0 0 U 3 0 0 0 0 0 0 0 0 U 3             , B =             0 0 0 E 3 (t) E 2 (t) -E 1 (t) 0 0             . ( 68
)
Here θ K and θ G are such that ∂ X3 θ = θ K ∂ X3 K + θ G ∂ X3 G, θ being defined by [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]. The gyrofluid model is hyperbolic: A(W) is diagonalizable and presents the eigenvalue λ 0 = U 3 of multiplicity 6 and two eigenvalues λ -= U 3 -c and λ + = U 3 + c with c = θ(1 + θ K ) of multiplicity 1. The associated eigenvectors for λ 0 are 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, θ G , -θ K ), r 6 = (-nθ K , 0, 0, 0, 0, 0, θ, 0). ( 69)

The other eigenvectors are given by r -= (n, 0, 0, -c, 0, 0, θ, 0) for λ -, r + = (n, 0, 0, +c, 0, 0, θ, 0) for λ + .

Remark that similarly to the classical Euler equations, the gyrofluid model presents two genuinely nonlinear fields associated to λ ± = U 3 ± c, as well as a linearly degenerated field associated to λ 0 = U 3 . To prove Proposition (2), we need the following Lemma.

Lemma 1. We have the relations below

∂ t n U 2 2 + ∂ X3 nU 3 U 2 2 = 0, (70)
∂ t n Y 2 2 + ∂ X3 nU 3 Y 2 2 = nY • ⊥ E, (71)
∂ t n U 2 2 + ∂ X3 nU 3 U 2 2 = U 2 2 (∂ t n + ∂ X3 nU 3 ) + n(U • (∂ t U + U 3 ∂ t 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

∂ t n + U 3 ∂ X3 n + n∂ X3 U 3 = 0 and ∂ t U + U 3 ∂ X3 U = 0,
whereas (65c) gives

∂ t U 3 + U 3 ∂ X3 U 3 + ∂ X3 θ + θ n ∂ X3 n = E 3 ,
where θ is a function of K and G which are independent of the other variables (see [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] and ( 19)). Then, there exist θ K (K, G) and θ G (K, G) such that ∂ X3 θ = θ K ∂ X3 K + θ G ∂ X3 G. Consequently, (65c) becomes

∂ t U 3 + U 3 ∂ X3 U 3 + θ K ∂ X3 K + θ G ∂ X3 G + θ n ∂ X3 n = E 3 . (73)
From (65d), we obtain

∂ t Y + U 3 ∂ X3 Y = ⊥ E.
We use (72) of Lemma 1 together with (65e) to get

∂ t nK + ∂ X3 nU 3 K + nθ∂ X3 U 3 = 0, or again ∂ t K + U 3 ∂ X3 K + θ∂ X3 U 3 = 0. ( 74
)
We use (70) and (71) of Lemma 1 to recast (65f) into

∂ t G + U 3 ∂ X3 G = 0.
We then deduce that (65) can be written in the quasilinear form (67). Moreover, we can check that θ K > 0 which implies that c = θ(1 + θ 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 -λId) = (U 3 -λ) 5 det

  U 3 -λ n 0 θ/n U 3 -λ θ K 0 θ U 3 -λ   .
It follows that det(A -λId) = (U 3 -λ) 6 (λ 2 -2U 3 λ + U 2 3 -c 2 ) = (U 3 -λ) 6 (λ -(U 3 -c))(λ -(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    

U 3 n 0 0 θ/n U 3 θ K θ G 0 θ U 3 0 0 0 0 U 3     .
This concludes the proof.

1 τ

 1 t) dt, and perform different quadratures (Lawson and exponential integrators type strategies) on the two last right hand side terms as follows Ŵn+1 p,l = e -( il ε +

Ŵn+1

  )-(48), we obtain when ε → 0

  L ∞ -Error for ε = 1 and different τ .

Figure 1 :

 1 Figure 1: Case M(w, z).

Figure 2 :

 2 Figure 2: Case M(z): L ∞ -Error for ε = 1 and different τ .

Figure 3 :

 3 Figure 3: Time history of the M(z(t)) and of ω(t) for different values of ε and τ .

Figure 4 :

 4 Figure 4: Weights trajectories for ∆t bigger than τ .

Figure 5 :

 5 Figure 5: Particle density relaxing to a gyromaxwellian.

3 Figure 6 :

 36 Figure 6: Evolution of the three errors e n i , i = 1, 2, 3 for scheme A for two different τ .

Figure 7 :

 7 Figure 7: Evolution of the three errors e n i , i = 1, 2, 3 for scheme B for two different τ .

nU 3 E 3 -U 3 ∂

 333 X3 nθ. (72) Proof. Starting from (65b) and use the continuity equation (65a) enables to get ∂ t U + U 3 ∂ X3 U = 0. Consequently, developing (70) leads to

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.


Aurora Rossi
email: aurora.rossi@inria.fr

Samuel Deslauriers-Gauthier

Emanuele Natale

Temporal Hyperbolic Graphs as Null Models for Brain Dynamics

Introduction

Null models are crucial for determining the degree of significance when testing hypotheses about brain dynamics modeled as a temporal complex network [START_REF] Lurie | Questions and controversies in the study of timevarying functional connectivity in resting fMRI[END_REF][START_REF] Váša | Null models in network neuroscience[END_REF]. The comparison between the hypothesis being tested on empirical data and on the null model enables us to assess the extent to which an apparently remarkable feature of the former can be attributed to randomness.

In this work, we consider networks generated from resting-state functional magnetic resonance imaging (fMRI) signals from the Human Connectome Project [START_REF] Wu-Minn | Subjects Data Release Reference Manual[END_REF] using standard processing techniques [START_REF] Preti | The dynamic functional connectome: State-of-the-art and perspectives[END_REF]. The network nodes correspond to brain regions while the edges represent the presence of a correlation greater than a certain value between the signals of two regions. Over time, edges may appear and disappear, indicating that the associated correlation may jump above and below a certain value.

While null models for static networks have been studied extensively [START_REF] Váša | Null models in network neuroscience[END_REF][START_REF] Fornito | Fundamentals of brain network analysis[END_REF][START_REF] Bassett | Small-World Brain Networks Revisited[END_REF][START_REF] Liao | Small-world human brain networks: Perspectives and challenges[END_REF], there is a lack of attention to temporal networks. Therefore, our investigation focuses on the study of temporal null models. In doing so, we focus on the ability of the null models to reproduce the temporal small-worldness present in the empirical data, a property associated with the efficient exchange of information at local and global scales over time.

Specifically, we compare two models taken from the literature [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF], which we will refer to as the Random Temporal Edge Graph Model and Random Temporal Permuted Times Graph Model, and three models that have not been previously considered in the context of brain dynamics which are the Random Temporal Square Graph Model, the Random Temporal Torus Graph Model [START_REF] Penrose | Random Geometric Graphs[END_REF] and the Random Temporal Hyperbolic Graph Model.

The latter is a temporal version of the Random Hyperbolic Graph Model which has gained prominence in the complex network research community for its ability to capture some important properties of real-world networks, such as small-worldness and a high-tail degree distribution [START_REF] Hartle | Dynamic Hidden-Variable Network Models[END_REF]. In addition, a recent paper shows that brain networks are best represented in the hyperbolic disk, which is a projection of hyperbolic space, showing relevance to its anatomical counterpart [START_REF] Whi | Hyperbolic disc embedding of functional human brain connectomes using resting-state fMRI[END_REF]. Our analysis indicates that it is also a particularly suitable null model for reproducing the temporal small-worldness observed in the empirical data.

Models

We consider five random temporal graph models. The Random Temporal Permuted Times Graph (RTPT) Model is created by randomizing the timing of each contact in the real data network [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF]. Similarly, the Random Temporal Edge Graph (RTE) Model is a permutation of the real data network where, in each snapshot, one end of an edge is reassigned to a different node [START_REF] Sizemore | Dynamic graph metrics: Tutorial, toolbox, and tale[END_REF].

The Random Temporal Square Graph (RTS) Model and the Random Temporal Torus Graph (RTT) Model are created by randomly placing nodes within a square and a torus, respectively and then connecting them if they are within a certain distance [START_REF] Penrose | Random Geometric Graphs[END_REF]. To introduce a temporal dependency to the models, the node positions are updated at each time step by applying a displacement vector with a uniformly chosen direction in [0, 2π) and length in (0, v), where v is the speed parameter of the models, and recomputing the edges accordingly.

The Random Temporal Hyperbolic Graph (RTH) Model is a temporal variant of Krioukov's Random Hyperbolic Graph (RH) model [START_REF] Hartle | Dynamic Hidden-Variable Network Models[END_REF]. Geometrically, nodes are quasiuniformly distributed within a disk with a radius of R, which is centered on the upper half of a hyperboloid with a specified negative curvature of K = -ζ 2 . Formally, writing the node positions in polar coordinates (r, θ ), where r ∈ [0, R] and θ ∈ [0, 2π], we have that the angular density is ρ(θ ) = 1 2π and the radial coordinate density is ρ(r) = α sinh(αr) cosh(αR)-1 . The α parameter controls the spread of the point positions. If α = ζ , the points follow a uniform distribution. If α > ζ , the points are more likely to be near the border of the disk, otherwise, they are more likely to be near the center. Two nodes, x and y, are connected if they are within hyperbolic distance d H < R, where, considering the polar coordinate (r x , θ x ) of x and (r y , θ y ) of y, d H is defined as To obtain the temporal version, we update the point position for as many time steps as there are in the real data networks we want to compare to. The point position update is chosen to ensure that the marginal distribution at each time step is the same as the initial distribution. Specifically, the polar coordinate is updated by adding to the previous value a number chosen according to a uniform distribution over (0, v), and then computing the modulo 2π. For the radial coordinate, we add to the previous value a number chosen according to a uniform distribution over (-v, v), and we reflect the result in the interval Temporal Small-Worldness The first static definition of small-worldness of a graph G was introduced by Watts and Strogatz [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF] as the combination of a high clustering coefficient C G (as in lattice graphs) and a short average path length L G (as in Erdős-Rènyi graphs). The previous measures are combined into the ratio C G L G . If the ratio is greater than 1, it indicates a significant small-world property of the network. To apply this to temporal networks, it is necessary to modify the previous measures to take into account the time evolution. We thus define the temporal small-worldness of a temporal graph

G T S G T = TC G T T L G T
where TC is the average of the global clustering coefficient over time and T L is defined as the average of the shortest path length between all pairs of nodes over time.

Optimization of the parameters The model parameters are determined using an optimization method that minimizes the area between the temporal small-worldness curves of the empirical data and the random temporal model over all different average degrees (see Fig. 1). The curves are created by interpolating the temporal small-worldness values at different average degrees. The average degrees differ based on the correlation threshold variation in the cases of empirical networks and the RTPT and RTE models, as well as the variation of the connection radius in the RTS, RTT, and RTH models.

Results

In Fig. 1, the blue lines are the values of the empirical network of 1042 subjects, created with the Schaefer atlas of 302 regions [START_REF] Schaefer | Local-global parcellation of the human cerebral cortex from intrinsic functional connectivity mri[END_REF]. For all the models, not depending on the real data, the illustrated lines are the median of the temporal small-worldness values over 10 random model realizations, instead, in the case of the real data, the RTE and RTPT models, the lines correspond to the median over 1042 subjects. The shadows correspond to the values between the first and third quartiles. Note that for RTT, RTS and RTH graph models the shadows are hardly visible given the minimal variance.

Looking at the random models, the median of the RTPT Graph Model has higher values with respect to those of the RTE Graph Model. We can observe that the RTT and RTS models perform well on high-average degree regimes. However, the RTH graph model not only exhibits the closest temporal small-worldness values but also follows the same trend as the empirical data.

d H (x, y) = 1 ζ

 1 acosh(cosh ζ r x cosh ζ r ysinh ζ r x sinh ζ r y cos(π -|π -|θ xθ y ||)).

[0, 1 ]

 1 to keep the distribution uniform in [0, 1]. As for the curvature parameter ζ , we set it to ζ = 1, since (Lemma 1.1 of [3]) if two RH graphs have the same ratio ζ α and the other parameters are equal, then they produce the same distribution on graphs. The parameters of the model are v and α.

Fig. 1 .

 1 Fig. 1.

Acknowledgments

This work has been supported by the French government, through the UCA DS4H Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-17-EURE-0004. The authors are grateful to the OPAL infrastructure from Université Côte d'Azur for providing resources and support.


Léa Roussin
email: lea.roussin@inrae.frwww.inrae.fr

Ines Haddam

Naika Prince

Lucía Peralta Marzal

Christine Heberden

Elise Maximin

Magali Monnoye

Catherine Phillipe

Paula Perez Pardo

Aletta Kraneveld

Sylvie Rabot

Laurent Naudon

Effect of gut microbiota from children with autism spectrum disorder on behavior and ASD-related biological markers in germ-free mice

Results: 16S rRNA sequencing analysis showed an increase in diversity in the AG group compared to S-AG according to three different indexes (Richness observed p= 0.0007; Shannon p=0.0011; Inversed Simpson p=0.013) (Fig. 1), accompanied by an increase in the relative abundance of the Desulfobacterota and Proteobacteria phyla (p=0.0459 and 0.0219 respectively) (Fig2). No difference of diversity was found between the S-A and A groups, however, the relative abundance of the Actinobacteria phylum was increased in the A group (p=0.0386) (Fig2). Many differences were also observed at lower phylogenetic levels, with more than 100 OTUs differing between each groups (not shown).

In parallel, analysis of the relative proportions of SCFAs in caecal content of the same animals (collected at sacrifice), showed a higher % of acetate in the AG group compared to S-AG (p=0.0156) as well as higher butyrate in both A and AG groups compared to S-A (p=0.0010) and S-AG (p=0.0225) respectively. A decrease of propionate was also observed in both A (p=0.0008) and AG (p=0.0004) groups compared to their respective sibling groups. Finally, branched and long chain SCFAs were decreased in A compared to S-A, (p=0.0066) but increased in AG compared to S-AG (p=0.253) (Fig. 3) Results: Recognition index (RI) in mice from the AG group is significantly lower than the one of the S-AG group (p=0.0030), indicating a decreased short term spatial memory in this group. In the A group however, RI was significantly higher than the one of the S-A group (p=0.261). In addition, only the A and S-AG groups showed a RI significantly different from the theoretical value of 50%, indicating a preference for the displaced objects (A: p=0.009; S-AG p=0.0004). These results indicate that the microbiota from children with ASD and GI symptoms negatively impacts spatial memory when implanted in originally GF C57BL/6 mice whereas microbiota from children with ASD and no GI symptoms seems to positively impact it. Although there is a strong genetic heritability to ASD, it is now accepted that environmental factors can play a role in its onset. As the prevalence of gastrointestinal (GI) symptoms is four-times higher in ASD patients, the potential implication of the gut microbiota in this disorder is being increasingly studied. A disturbed microbiota composition has been demonstrated in ASD patients, accompanied by altered production of bacteria metabolites. This project aims to determine if the microbiota can influence behavior when transplanted to germ-free mice. Our hypothesis is that a pool of microbiota from patients with ASD, especially those with GI symptoms, will worsen behavioral phenotypes related to ASD in those mice as well as impact various ASD related biological markers, in comparison to GF mice transplanted with microbiota for the neurotypical siblings of those children. To ensure that the mice were exclusively exposed to donor microbiota, they were kept in isolators for the duration of the experiments. Experiments have been carried out on both BALB/c and C57BL/6 mice, but some analyses are still ongoing.

Sas

  Fig 1: α-diversity measures of the microbiota (mice feces 6 weeks after FMT)

Fig 3 :Fig 5 :

 35 Fig 3: Proportion of SCFAs in mice caecal content (analysed by gaz chromatograpy)

  (ASD) is a neurodevelopmental disorder affecting 1 in 160 people in the world.

Acknowlegments:

Thanks to the Migale platform for their script, tools, and advice for 16S analysis. Thanks to Patricia Alves Soares and Johana Jarkulischová for their work on the immunohistochemistry. This work is part of the GEMMA project funded by the European commission (Horizon 2020 program, grant agreement ID 825033) Images from biorender.com

S-A A S-AG AG Conclusion

The microbiota that was transferred to the BALB/c mice has a distinct diversity, composition and fermentation activity in our four groups, which is consistent with literature for ASD patients. For C57BL/6 mice, the analysis is still ongoing.

The effect of the FMT on behavior was limited. In BALB/c mice we did not observe any difference between our four groups. In C57BL/6 mice, we saw an effect of the microbiota on short term spatial memory, however, it is hard to interpret this result are the effects are opposite between A and AG groups. It may be noted that two previous papers (Sharon et al. Cell, 2019;and Xiao et al. ASM, 2021) with similar study designs both found an impact of microbiota from patients with ASD on anxiety, social and repetitive behaviors in GF mice.

However, the different microbiota did have an impact of biochemical factors, as the number of serotoninergic or serotonin positive neurons is decreased in the brain and the gut in both "ASD" groups in BALB/c mice. This is interesting as the serotonin pathway has been shown to be dysregulated, both in ASD patients and various animal models of ASD (Muller et al. Neuroscience, 2016). Those analyses are still ongoing for C57BL/6 mice, and further analyses on gut inflammation and permeability are still ongoing for both strains, which will hopefully bring more information on the effects of the FMT at a biochemical level. Overall, our results show that the impact of the gut microbiota on behavioral manifestations of ASD is still unclear, and more research is needed to further this hypothesis.


Autism spectrum disorder embracingasd.com tacanow.org

• High prevalence of gastrointestinal (GI) disorders in ASD patients

• Genetic heritability is reported to be around 50% also a strong environmental impact.

Role of microbiota-gut-brain axis in ASD ?

• Several studies observe differences in the composition of the microbiota between people with ASD and neurotypical people. The presence of GI symptoms seems to also have an impact.

• Little inter-study consensus on the bacterial genera or species modulated, probably due to differences in geography, diet, and choice of control groups (siblings or not).

•A few studies show improvment of behavior and/or GI symptoms after microbiota modulation (FMT, probiotics) • Inter-group differences in caecal SCFA profiles

••→

  in ASD patients or mice models of ASD Less serotonin positive cells in the ileum of A and AG mice in BALBc and a trend for group Reduced number of serotonergic neurons in raphe nuclei in groups A and AG only Lower recognition index in S-A than A mice •No difference in anxiety, social behavior and stereotyped behavior

BALB/cC57BL/6

Conclusion:

• A more in depth statistical analysis is planned •Some analysis are still ongoing (quantification of gene expression in the brain)

Microbiota