Mihai Bostan
email: mbostan@univ-fcomte.fr

Claudia Negulescu
email: claudia.negulescu@cmi.univ-mrs.fr

Mathematical models for strongly magnetized plasmas with mass disparate particles

Keywords: Vlasov equation, Multi-scale analysis, Average operator AMS classification: 35Q75, 78A35, 82D10

The controlled fusion is achieved by magnetic confinement : the plasma is confined into toroidal devices called tokamaks, under the action of strong magnetic fields. The particle motion reduces to advection along the magnetic lines combined to rotation around the magnetic lines. The rotation around the magnetic lines is much faster than the parallel motion and efficient numerical resolution requires homogenization procedures. Moreover the rotation period, being proportional to the particle mass, introduces very different time scales in the case when the plasma contains disparate particles; the electrons turn much faster than the ions, the ratio between their cyclotronic periods being the mass ratio of the electrons with respect to the ions. The subject matter of this paper concerns the mathematical study of such plasmas, under the action of strong magnetic fields.

In particular, we are interested in the limit models when the small parameter, representing the mass ratio as well as the fast cyclotronic motion, tends to zero.

Introduction

Many research programs in plasma physics are devoted to magnetic confinement. It concerns the dynamics of a population of charged particles under the action of strong magnetic fields, let say B ε (x) depending on some parameter ε > 0. For instance, thermonuclear fusion (and thus energy) is produced in a tokamak, which is a toroidal plasma confining device, the ionized gaz (plasma) being confined by a strong magnetic field. Indeed, the radius of the circular motion of the charged particles around the magnetic lines (which is called the Larmor radius ρ L ) is proportional with the inverse of the magnetic field, i.e. ρ L = mv/|qB ε |. Here m is the particle mass, q is the particle charge and v is the velocity in the plane perpendicular to the magnetic field lines. Therefore, when the magnetic field becomes large, the Larmor radius vanishes and thus, at the lowest order, the particles remain confined around the magnetic lines, which are supposed to enclose a bounded volume (the tokamak). But strong magnetic fields introduce small time scales, since the rotation period of particles (called cyclotronic period T c ) is proportional to the inverse of the magnetic field, i.e. T c = 2πm/|qB ε |.

Clearly, the efficient numerical resolution of such models requires multiple scale analysis or homogenization techniques. Notice also that in the case of a gaz, consisting of distinct particles (let us say ions/electrons), the cyclotronic motion introduces several small time scales, for example when the particle masses are disparate.

Using the kinetic description and neglecting the collisions we are led to the Vlasov equations

∂ t f ε i + p m i • ∇ x f ε i + e E(t, x) + p m i ∧ B ε (x) • ∇ p f ε i = 0, (t, x, p) ∈ R + × R 3 × R 3 (1)
∂ t f ε e + p m e • ∇ x f ε e -e E(t, x) + p m e ∧ B ε (x) • ∇ p f ε e = 0, (t, x, p) ∈ R + × R 3 × R 3 (2)
where f ε i (resp. f ε e ) is the distribution function of the ions (resp. electrons) in the position-momentum phase space (x, p) ∈ R 3 × R 3 , m i (resp. m e ) is the ion (resp. electron) mass and e (resp. -e) is the ion (resp. electron) charge. Remark that the system (1), [START_REF] Allen | Extended-geostrophic Hamiltonian models for rotating shallow water motion[END_REF] is not written in the usual position-velocity phase space, as we are assuming for the moment, that the electron and ion momenta are of the same order of magnitude.

We assume that B ε is a given stationary, divergence free magnetic field and that the electric field derives from a given electric potential E(t, x) = -∇ x φ(t, x). We suppose moreover that the electro-magnetic field is smooth. We are interested now in the behaviour of the system (1), (2) for a large magnetic field and also large mass ratio between ions and electrons m i m e = 1 µ >> 1.

Choosing T ci as a typical value of the set {T ci (x) : x ∈ R 3 }, we assume that the time observation is much larger than the typical ion cyclotronic period

T obs T ci = 1 ε >> 1 .
This leads to a magnetic field of the form

B ε (x) = B(x) ε b(x), div x (Bb) = 0
for some scalar positive function B(x) (given by eB(x)/m i = 2πT ci /(T obs T ci (x)), x ∈ R 3 ) and some field of unitary vectors b(x). Introducing moreover the rescaled ion cyclotronic frequency

ω ci (x) = eB(x) m i
we obtain from (1), [START_REF] Allen | Extended-geostrophic Hamiltonian models for rotating shallow water motion[END_REF] the Vlasov equations

∂ t f ε i + p m i • ∇ x f ε i + e E(t, x) + 1 ε ω ci (x)p ∧ b(x) • ∇ p f ε i = 0 (3)
∂ t f ε e + 1 µ p m i • ∇ x f ε e -e E(t, x) + 1 µε ω ci (x)p ∧ b(x) • ∇ p f ε e = 0 (4)
where ε and µ are small parameters relating the observation time and the typical ion/electron cyclotronic periods

T obs T ci = 1 ε , T ci T ce = m i m e = 1 µ , T obs T ce = 1 εµ .
Certainly, the equations ( 3), (4) can be written in dimensionless form, by introducing the unknowns and variables f ε i/e ( t, x, p) = 1 X 3 obs P 3 obs f ε i/e (t, x, p), Ẽ( t, x) = P obs eT obs E(t, x) t = T obs t, x = X obs x, p = P obs p, P obs m i = X obs T obs .

In this case we obtain

∂ t f ε i + p • ∇ x f ε i + Ẽ( t, x) + 2π ε T ci T ci (x) p ∧ b(x) • ∇ p f ε i = 0 (5) ∂ t f ε e + p µ • ∇ x f ε e -Ẽ( t, x) + 2π µε T ci T ci (x) p ∧ b(x) • ∇ p f ε e = 0. (6)
Nevertheless, in order to obtain the standard quantities, like magnetic momentum, electric cross field drift, magnetic curvature/gradient drift expressed in physical units, we prefer to work with the equations ( 3), ( 4), instead of ( 5), [START_REF] Bostan | Gyrokinetic Vlasov equation in three dimensional setting[END_REF].

The ordering between these time scales is

T obs >> T ci >> T ce .
We analyze here the particular regime characterized by T ci = √ T obs T ce i.e., when ε = µ. Generally the ratio between the observation time and the ion cyclotronic period is chosen as a function of the ratio between the typical ion Larmor radius ρ Li and the tokamak small radius a, T obs /T ci = (a/ρ Li ) 2 . Since T ci /T ce = m i /m e , then the regime

T ci = √
T obs T ce leads to a/ρ Li = m i /m e = √ 1836 (which corresponds to tokamaks like TFTR Princeton). Having dropped the index i for the ion mass and the rescaled cyclotronic frequency, the Vlasov equations to be studied become

∂ t f ε i + p m • ∇ x f ε i + e E(t, x) + 1 ε ω c (x)p ∧ b(x) • ∇ p f ε i = 0 ( 7
)
∂ t f ε e + 1 ε p m • ∇ x f ε e -e E(t, x) + 1 ε 2 ω c (x)p ∧ b(x) • ∇ p f ε e = 0. (8)
The aim of this work is to investigate the asymptotics of ( 7), [START_REF] Boozer | Physics of magnetically confined plasmas[END_REF] when ε ց 0 and to emphasize the specific behaviour of the ions and electrons under the action of a strong magnetic field. Notice that doing this, we assume that the electrons and ions posses momenta of the same order of magnitude. Later on (see Section 4.3), we shall investigate the asymptotics under the hypothesis that electrons and ions have similar velocities. The two limit models are rather different and correspond to regions of the tokamak with disparate temperatures between ions and electrons, as illustrated by the formula

θ i θ e = m i m e v th i v th e 2 = m e m i p th i p th e 2
where θ i /θ e , v th i /v th e and p th i /p th e are the temperatures, thermal velocities and thermal momenta of ions/electrons. As usual in multiple scale analysis, the main idea is to separate scales, i.e. to distinguish between fast and slow motion. The effective dynamics is obtained by averaging over the small time scale, here the cyclotronic period. As observed in [START_REF] Bostan | Transport equations with singular coefficients. Application to the gyrokinetic models in plasma physics[END_REF] this approach can be interpreted from the ergodic point of view: it reduces to mean ergodic theorem, which allows us to construct an average operator, associated to the smallest time scale.

The goal of this paper is how to generalize this method when two different small time scales appear in the model, as for example in the electron Vlasov equation [START_REF] Boozer | Physics of magnetically confined plasmas[END_REF]. The main idea is to perform double average, one for each small time scale. More exactly we start by averaging with respect to the smallest time scale ε 2 . A new Vlasov equation is obtained, with only one small time scale left. Finally a second average is performed, in order to remove the fluctuations evolving on the time scale ε. Up to our knowledge, this method which combines successive average operators in order to handle several small scales is new.

If the computations are completely explicit for ions and this for general magnetic shapes, things are more complex for electrons. We obtain explicit formula at least in some particular cases (cylindrical geometry), as

b(x) = (x 2 , -x 1 , 1) (x 2 1 + x 2 2 + 1) 1/2 , x ∈ R 3 . (9)
More generally, the arguments presented in this paper allow to treat many other models, not only the case of strongly magnetized plasmas with disparate particle masses. The method can be adapted straightforwardly to any linear transport equation involving multiple scales ε, ε 2 , ..., ε p with p ∈ N ⋆ , but the explicit derivation of the limit model may become very complex since, in general, it requires p averaging processes.

The goal of this paper is to distinguish the dynamics of ions and electrons by taking into account their mass ratio. Two different ion/electron limit models will be obtained, depending on the starting assumption of similar ion/electron momenta or velocities.

The relative mass constraint between two particle species has been addressed in previous works devoted to kinetic theory (Boltzmann equation, Fokker-Planck equation) [START_REF] Degond | Transport coefficients of plasmas and disparate mass binary gases[END_REF], [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF], [START_REF] Charles | Small mass ratio limit of Boltzmann equations in the context of the study of evolution of dust particles in a rarefied atmosphere[END_REF] but not under the hypothesis of strong magnetic field. The new contribution of the present work is to provide a rigorous mathematical analysis which describes the magnetic confinement of several species of charged particles and explain the specific behaviour when keeping trace of their relative mass.

For the analysis of the Vlasov or Vlasov-Poisson equations with a large external magnetic we mention [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF], [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF], [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF], [START_REF] Frénod | The finite Larmor radius approximation[END_REF]. The numerical approximation of the gyrokinetic models has been performed in [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF] using semi-Lagrangian schemes.

The nonlinear gyrokinetic theory of the Vlasov-Maxwell equations can be carried out by appealing to Lagrangian and Hamiltonian methods [START_REF] Boozer | Physics of magnetically confined plasmas[END_REF], [START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF], [START_REF] Littlejohn | A guiding center Hamiltonian : A new approach[END_REF], [START_REF] Littlejohn | Hamiltonian formulation of guiding center motion[END_REF]. It is also possible to follow the general method of multiple time scale or averaging perturbation developped in [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations[END_REF]. For a unified treatment of the main physical ideas and theoretical methods that have emerged on magnetic plasma confinement we refer to [START_REF] Hazeltine | Plasma confinement[END_REF], [START_REF] Rax | Physique des plasmas, Cours et applications[END_REF].

We also mention that the drift approximation of strongly magnetized plasmas is analogous to the geostrophic flow in the theory of a shallow rotating fluid [START_REF] Abarbanel | Hamiltonian description of almost geostrophic flow[END_REF], [START_REF] Allen | Extended-geostrophic Hamiltonian models for rotating shallow water motion[END_REF], [START_REF] Chemin | Mathematical geophysics: An introduction to rotating fluids and the Navier-Stokes equations[END_REF], [START_REF] Salmon | New equations for nearly geostrophic flow[END_REF], [START_REF] Vanneste | Dirac-bracket approach to nearly geostrophic Hamiltonian balanced models[END_REF].

Our paper is organized as follows. Section 2 presents briefly the main ideas of this work, as well as the main results. In Section 3 we introduce the first average operator and list its mathematical properties : orthogonal decomposition of L 2 functions into zero average functions and invariant functions along the characteristic flow, Poincaré inequality, etc. The ion limit model follows immediately by averaging along the characteristic flow corresponding to the dominant transport operator in [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF]. In the first part of Section 4 we introduce the second average operator, since the analysis of the electron limit model ( 8) requires double averaging. The second part of Section 4 is devoted to the asymptotics of the electron model (8) (momentum units of the same order for both ions and electrons). We investigate magnetic shapes whose lines are winding on cylindrical surfaces. In the last part of Section 4 we perform a similar analysis when assuming that the velocity units are of the same order for both ions and electrons.

Presentation of the models and main results

The dynamics in [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] is dominated by the transport operator

1 ε ω c (x)(p ∧ b(x))
• ∇ p and leads to the guiding-center approximation cf. [START_REF] Bostan | Gyrokinetic Vlasov equation in three dimensional setting[END_REF]. A formal derivation follows by using a standard asymptotic expansion like

f ε i = f i + εf 1 i + ε 2 f 2 i + ... (10)
Plugging the above Ansatz in [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] and denoting by T the operator ω c (p∧b)•∇ p yields at the lowest order the divergence constraint T f i = 0 and to the next order the evolution equation for f i

∂ t f i + p m • ∇ x f i + e E(t, x) • ∇ p f i + T f 1 i = 0. ( 11
)
We need to close [START_REF] Chemin | Mathematical geophysics: An introduction to rotating fluids and the Navier-Stokes equations[END_REF] with respect to the first order fluctuation density f 1 i . Motivated by the fact that the leading order term f i belongs to the kernel of T , we project [START_REF] Chemin | Mathematical geophysics: An introduction to rotating fluids and the Navier-Stokes equations[END_REF] on ker T . Since the range of T is orthogonal to its kernel, it is easily seen that in this way we can eliminate the unknown f 1 i from ( 11) and obtain a transport equation which permits to compute f i . For the explicit computation of the advection field of this Vlasov like limit equation, it is worth observing that the orthogonal projection on ker T is equivalent to averaging along the characteristic flow associated to T . The rigorous construction of the average operator (sometimes called by physicists the gyro-average operator in the context of gyrokinetic models) essentially relies on ergodic theory i.e., von Neumann's ergodic theorem ([22] pp. 57). Employing this method yields the Vlasov equation ( 16) (see Proposition 2.1) for the ion dominant term f i , cf. [START_REF] Bostan | Gyrokinetic Vlasov equation in three dimensional setting[END_REF]. If the derivation of the limit model ( 16) is now well understood, the behaviour of (8) when ε ց 0 is not obvious. Assuming that the electron density f ε e for small ε, behaves like

f ε e = f e + εf 1 e + ε 2 f 2 e + ... (12)
leads to the following equations corresponding to the orders ε -2 , ε -1 , ε 0 , ...

T f e = 0 ( 13
) p m • ∇ x f e -T f 1 e = 0 ( 14
)
∂ t f e -e E • ∇ p f e + p m • ∇ x f 1 e -T f 2 e = 0. ( 15
)
The key point is how to close [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF] with respect to the first and second order fluctuation terms f 1 e , f 2 e . Certainly the constraints [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF] and [START_REF] Frénod | The finite Larmor radius approximation[END_REF] have to be taken into account. We intend to perform our analysis using average operators, as it was done for the ion dynamics. The problem is more complex since (8) combines two different scales : ε and ε 2 . We will see that the limit model can be obtained by similar techniques, involving double averaging. In the particular case of a magnetic shape whose lines are winding on cylindrical surfaces [START_REF] Hazeltine | Plasma confinement[END_REF] an explicit Vlasov equation is derived for the dominant electron density. The ion/electron limit models obtained under the assumption of similar ion/electron momenta are summarized in Proposition 2.1 Let us assume that the electro-magnetic field is smooth, the magnetic field B ε being divergence free and of the form

B ε (x) = B(x) ε b(x) , |b(x)| = 1 , x ∈ R 3
whereas the electric field is given via a potential, as E(t, x) = -∇ x φ(t, x). Moreover, let us assume that inf x∈R 3 B(x) > 0.

i) Ion limit model

Then the limit ion density f i = lim εց0 f ε i , with f ε i solving (7), satisfies

∂ t f i + b(x) ⊗ b(x) p m • ∇ x f i + eb(x) ⊗ b(x)E + ω i (x, p) ⊥ p • ∇ p f i = 0 ( 16
)
where for any (x, p) with p∧b(x) = 0 the symbol ⊥ p stands for the orthogonal momentum to p, contained in the plane determined by b(x) and p, and such that its coordinate along b(x) is positive, that means

⊥ p = |p ∧ b(x)| b(x) -(p • b(x)) b(x) ∧ (p ∧ b(x)) |p ∧ b(x)| (17)
and the frequency ω i (x, p) is given by

ω i (x, p) = |p ∧ b(x)| 2m div x b - (p • b(x)) m ∂b ∂x b(x) • p |p ∧ b(x)| , p ∧ b(x) = 0.
ii) Electron limit model

We shall assume in this case the particular magnetic shape

b(x) = (x 2 , -x 1 , 1) (x 2 1 + x 2 2 + 1) 1/2 , x ∈ R 3 . ( 18
)
Then the limit electron density f e = lim εց0 f ε e , with f ε e solving (8), satisfies

∂ t f e + (v GD + v CD ) • ∇ x f e + ω e (x, p) ⊥ p • ∇ p f e = 0
where the gradient drift resp. curvature drift velocities are defined as

v GD = |p ∧ b| 2 2m 2 (-ω c ) b ∧ ∇ x B B , v CD = (p • b) 2 m 2 (-ω c ) (b ∧ ∂ x b b)
and the electron frequency ω e is given by

ω e (x, p) = -(v GD + v CD ) • t ∂ x b p |p ∧ b| + |p ∧ b|(p • b) 2m 2 ω c ∇ x B B • (b ∧ ∂ x b b).
Notice that the particular case of Proposition 2.1 ii) still captures the main drifts, the magnetic curvature/gradient drifts, as predicted for general magnetic shapes.

Another interesting asymptotic case is that of typical ion/electron velocities of the same order of magnitude. Denoting by F ε i (resp. F ε e ) the distribution function of the ions (resp. electrons) in the position-velocity phase space (x, v) ∈ R 3 × R 3 , the starting Vlasov equations are now (see (1), (2) for comparison)

∂ t F ε i + v • ∇ x F ε i + e m i (E(t, x) + v ∧ B ε (x)) • ∇ v F ε i = 0, (t, x, v) ∈ R + × R 3 × R 3 (19) ∂ t F ε e + v • ∇ x F ε e - e m e (E(t, x) + v ∧ B ε (x)) • ∇ v F ε e = 0, (t, x, v) ∈ R + × R 3 × R 3 . ( 20
)
The same ordering as previously T obs T ci = T ci Tce = m i me = 1 ε >> 1, leads to the models

∂ t F ε i + v • ∇ x F ε i + eE(t, x) m + 1 ε ω c (x) v ∧ b(x) • ∇ v F ε i = 0, (t, x, v) ∈ R + × R 3 × R 3 (21)
∂ t F ε e + v • ∇ x F ε e - 1 ε eE(t, x) m + 1 ε 2 ω c (x) v ∧ b(x) • ∇ v F ε e = 0, (t, x, v) ∈ R + × R 3 × R 3 (22)
where m = m i , ω c = ω ci . Clearly, the limit model for the ion density is similar to [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF].

When neglecting the parallel component of the electric field (i.e., E • b = 0) it happens that at the lowest order the electron density satisfies the same model as the ion density.

But the ions/electrons behave differently when first order corrections are taken into account. In particular the ions/electrons deviate differently from the magnetic lines.

With the notations

F ε i = F i + εF 1 i + ε 2 F 2 i + ..., F ε e = F e + εF 1 e + ε 2 F 2 e + ...
under the assumption of similar ion/electron velocities we prove Proposition 2.2 Assume that the electric and (rescaled) magnetic field are given and smooth, such that E(t, x) • b(x) = 0.

0) The zeroth order ion/electron densities F i resp. F e are solutions of the same limit model (similar to (16))

∂ t F i/e + b(x) ⊗ b(x)v • ∇ x F i/e + Ω i (x, p) ⊥ v • ∇ v F i/e = 0 ( 23
)
where

Ω i (x, p) = |v ∧ b(x)| 2 div x b(x) -(v • b(x)) ∂ x b b(x) • v |v ∧ b(x)| and ⊥ v = |v ∧ b(x)| b(x) -(v • b(x)) v -(v • b(x)) b(x) |v ∧ b(x)| .
i) The mean electron drift (up to second order corrections) is given by

R 3 (F e + εF 1 e )(v -(v • b) b) dv R 3 F e dv = ε E ∧ b B .
ii) The mean ion drift (up to second order corrections) is given by

R 3 (F i + εF 1 i )(v -(v • b) b) dv R 3 F i dv = ε E ∧ b B + V 2 ⊥ 2ω c b ∧ ∇ x B B + V 2 ω c - V 2 ⊥ 2ω c k(x)b(x) ∧ n(x) + b ∧ ∇ x R 3 µ(x, v)F i dv e R 3 F i dv where V ⊥ (t, x) = R 3 |v ∧ b(x)| 2 F i (t, x, v) dv R 3 F i (t, x, v) dv 1/2 , V (t, x) = R 3 (v • b(x)) 2 F i (t, x, v) dv R 3 F i (t, x, v) dv 1/2 and µ(x, v) = m|v ∧ b(x)| 2 /2B(x), k(x) = |∂ x b b(x)| is the curvature of the magnetic lines and n(x) = ∂ x b b(x)/k(x)
is the first normal to the magnetic lines.

3 First average operator and the ion limit model

The concern of this section shall be the introduction of the first average operator needed for the obtention of the ion/electron limit models, as ε → 0. We present in detail the ion model, in order to facilitate the understanding of Section 4. Recall that in this section, we suppose that the ions and electrons have momenta of the same order of magnitude.

First average operator

Our study is based on the construction of average operators, corresponding to characteristic flows preserving the Lebesgue measure. We work in the L 2 (R 3 × R 3 ) framework and we define the operator

T u = div p (ω c (x) u p ∧ b(x)) , u ∈ D(T ) D(T ) = {u(x, p) ∈ L 2 (R 3 × R 3 ) : div p (ω c (x) u p ∧ b(x)) ∈ L 2 (R 3 × R 3 )}.
The notation • stands for the standard norm of L 2 (R 3 × R 3 ). We denote by (X, P )(s; x, p) the characteristics associated to the vector field (0,

ω c (x)(p ∧ b(x))), that means dX ds = 0, dP ds = ω c (X(s)) P (s) ∧ b(X(s)), (X, P )(0) = (x, p). (24)
It is easily seen that x, |p ∧ b(x)|, (p • b(x)) are left invariant along the characteristic flow [START_REF] Vanneste | Dirac-bracket approach to nearly geostrophic Hamiltonian balanced models[END_REF]. Notice that each vector p ∈ R 3 can be decomposed into its parallel part p with respect to the magnetic field lines and its orthogonal part p ⊥ , like

p = p + p ⊥ , |p | 2 + |p ⊥ | 2 = |p| 2 , with p := (p • b(x)) b(x) = b(x) ⊗ b(x) p , p ⊥ := b(x) ∧ (p ∧ b(x)) = (I -b(x) ⊗ b(x))p ,
and where the symbol u ⊗ v, with u, v ∈ R 3 , stands for the matrix (u i v j ) 1≤i,j≤3 . The reader has to distinguish between the two different notations p ⊥ and ⊥ p given in [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF].

Straightforward computations yield the formulae X(s; x, p) = x and

P (s; x, p) = cos(ω c (x)s) p ⊥ + sin(ω c (x)s) p ⊥ ∧ b(x) + p .
The trajectories (X, P )(s; x, p) are T c (x) = 2π ωc(x) periodic for any initial condition (x, p) ∈ R 3 ×R 3 and therefore we introduce the average operator along these trajectories cf. [START_REF] Bostan | Transport equations with singular coefficients. Application to the gyrokinetic models in plasma physics[END_REF]

u (x, p) = 1 T c (x) Tc(x) 0 u(X(s; x, p), P (s; x, p)) ds = 1 2π S(x) u(x, |p ∧ b(x)| ω + (p • b(x)) b(x)) dω for any function u ∈ L 2 (R 3 × R 3 ), where S(x) = {ω ∈ S 2 : b(x) • ω = 0}. It is easily seen that Proposition 3.1
The average operator is linear and continuous. Moreover it coincides with the orthogonal projection on the kernel of T i.e.,

• :

L 2 (R 3 × R 3 ) → ker T and R 3 R 3 (u -u )ϕ dpdx = 0, ∀ ϕ ∈ ker T . Proof For any function u ∈ L 2 (R 3 × R 3 ) we have for a.a. x ∈ R 3 | u | 2 (x, p) ≤ 1 T c (x) Tc(x) 0 u 2 (x, P (s; x, p)) ds.
Taking into account that for any x ∈ R 3 the map p → P (s; x, p) is measure preserving one gets

R 3 R 3 u 2 (x, p) dpdx ≤ R 3 R 3 u 2 (x, p) dpdx saying that • ∈ L(L 2 (R 3 × R 3 ), L 2 (R 3 × R 3 )) and • L(L 2 (R 3 ×R 3 ),L 2 (R 3 ×R 3 )) ≤ 1.
It is well known that the kernel of T is given by the functions in L 2 invariant along the characteristics [START_REF] Vanneste | Dirac-bracket approach to nearly geostrophic Hamiltonian balanced models[END_REF]. Therefore we have

ker T = {u ∈ L 2 (R 3 × R 3 ) : ∃ v such that u(x, p) = v(x, |p ∧ b(x)|, (p • b(x)))} Notice that for any u ∈ L 2 (R 3 ×R 3 ) its average u depends only on x, |p∧b(x)|, (p•b(x)). Therefore u ∈ ker T . Pick a function ϕ ∈ ker T i.e., ∃ ψ : ϕ(x, p) = ψ(x, |p ∧ b(x)|, (p • b(x))) ∈ L 2 (R 3 × R 3 )
and let us compute

I = R 3 R 3 (u -u )ϕ dpdx. Using cylindrical coordinates along b(x) axis yields I = R 3 R R + ψ(x, r, z) S(x) u(x, r ω + z b(x)) dω -2π u rdrdzdx = 0 and therefore u = Proj ker T u for any u ∈ L 2 (R 3 × R 3 ). In particular u = u for any u ∈ ker T and • L(L 2 (R 3 ×R 3 ),L 2 (R 3 ×R 3 )) = 1.
The above result allows us to characterize the closure of the range of T . Indeed, since

• = Proj ker T and T ⋆ = -T we have ker • = (ker T ) ⊥ = (ker T ⋆ ) ⊥ = Range T .
Moreover we have the orthogonal decomposition of L 2 (R 3 ×R 3 ) into invariant functions along the characteristics [START_REF] Vanneste | Dirac-bracket approach to nearly geostrophic Hamiltonian balanced models[END_REF] and zero average functions i.e., L

2 (R 3 × R 3 ) = ker T ⊕ ⊥ ker • since u = u + (u -u ), R 3 R 3 (u -u ) u dpdx = 0.
If the magnetic field remains away from 0, the range of T is closed, leading to the equality Range T = ker • , which gives a solvability condition for T u = v. For the sake of the presentation we recall here the Poincaré inequality cf. where

ω 0 = e m inf x∈R 3 B(x) > 0.
Proof By the previous computations we know that Range T ⊂ ker

• . Assume now that u ∈ D(T ) ∩ ker • such that T u = 0. Since • = Proj ker T we have u = u = 0 saying that T | ker • is injective. Consider now v ∈ ker • and let us prove that there is u ∈ ker • ∩ D(T ) such that T u = v. For any α > 0 there is a unique u α ∈ D(T ) such that α u α + T u α = v. ( 25
)
Indeed it is easily seen that the solutions (u α ) α>0 are given by

u α (x, p) = R - e αs v(x, P (s; x, p)) ds, (x, p) ∈ R 3 × R 3 .
Applying the average operator to (25) yields u α = 0 for any α > 0. We are looking now for a bound of ( u α ) α>0 . We introduce the function V (s; x, p) = 0 s v(x, P (τ ; x, p)) dτ . Notice that for any fixed (x, p) the function s → V (s; x, p) is T c (x) periodic, because

v = 0 and thus V (s; x, •) L 2 (R 3 ) ≤ T c (x) v(x, •) L 2 (R 3 ) for any s ∈ R. Integrating by parts we obtain u α (x, p) = - R - e αs ∂ s V ds = R - αe αs V (s; x, p) ds implying that u α (x, •) L 2 (R 3 ) ≤ R - αe αs V (s; x, •) L 2 (R 3 ) ≤ T c (x) v(x, •) L 2 (R 3 ) ≤ T 0 v(x, •) L 2 (R 3 ) ,
where T 0 = 2π |ω 0 | . After integration with respect to x we obtain the uniform estimate u α ≤ T 0 v for any α > 0. Extracting a sequence (α n ) n such that lim n→+∞ α n = 0,

lim n→+∞ u αn = u weakly in L 2 (R 3 × R 3 ) we deduce easily that u ∈ D(T ), T u = v, u = 0, u ≤ T 0 v saying that T | ker • -1 is bounded linear operator and T | ker • -1 L(ker • ,ker • ) ≤ T 0 .

Ion limit model

Using the properties of the average operator • , we can easily derive the limit model ( 16) for the ion distribution, stated in Proposition 2.1 (i). Presenting a complete rigorous justification of the expansion [START_REF] Charles | Small mass ratio limit of Boltzmann equations in the context of the study of evolution of dust particles in a rarefied atmosphere[END_REF] is not one of the major priorities in this paper. The main objective is to provide a robust method for analyzing the asymptotics of linear transport equations like [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] and how to extend it to multi-scale problems like [START_REF] Boozer | Physics of magnetically confined plasmas[END_REF]. Nevertheless a rigorous weak convergence result is presented in Proposition 6.1.

We emphasize that the method we employ here has been studied in detail in [START_REF] Bostan | Transport equations with singular coefficients. Application to the gyrokinetic models in plasma physics[END_REF] (see also [START_REF] Bostan | The Vlasov-Poisson system with strong external magnetic field. Finite Larmor radius regime[END_REF]) for linear transport problems with even more general dominant advection fields, with characteristic flows not necessarily periodic. We refer to these papers for a complete mathematical analysis justifying rigorously the asymptotic behaviour.

Proof of Proposition 2.1 (i)

Let us recall that Ansatz (10) yields the evolution equation for the zeroth order distribution function f i

∂ t f i + p m • ∇ x f i + e E(t, x) • ∇ p f i + T f 1 i = 0 , (26)
under the constraint T f i (t) = 0. This constraint implies that there is a function such that we can write

g i = g i (t, x,
f i (t, x, p) = g i (t, x, |p ∧ b(x)|, (p • b(x))). (27)
Under the hypothesis in Proposition 3.2 the equality ker • = Range T holds true and therefore (26) is equivalent to

∂ t f i + p m • ∇ x f i + eE(t) • ∇ p f i = 0. ( 28
)
It remains to average the time, position and momentum derivatives of the dominant term f i . It is easily seen that the time derivative and the average operator are commuting since the characteristic system (24) is autonomous. Taking into account that

f i ∈ ker T we obtain ∂ t f i = ∂ t f i = ∂ t f i . (29)
For computing the averages of the space and momentum derivatives we apply the chain rule to (27) and we average only the derivatives of the invariants since the derivatives of g i depend only on time and the invariants and thus are constant along the characteristic flow [START_REF] Vanneste | Dirac-bracket approach to nearly geostrophic Hamiltonian balanced models[END_REF]. Assume for the moment that g i is smooth. By direct computations one gets for any (x, p) such that p ∧ b(x) = 0

p • ∇ x f i = p • ∇ x g i -∂ r g i (p • b(x)) |p ∧ b(x)| (∂ x b : p ⊗ p) + ∂ z g i (∂ x b : p ⊗ p) .
and

∇ p f i = ∂ r g i |p ∧ b(x)| (I -b(x) ⊗ b(x))p + ∂ z g i b(x).
Here the notation U : V stands for the contraction

3 i,j=1 u ij v ij of two matrices U = (u ij ), V = (v ij ) ∈ M 3×3 (R). Recall that p → (I -b(x) ⊗ b(x)
)p is the orthogonal projection and p → b(x) ⊗ b(x)p is the parallel projection with respect to the plane oriented by the magnetic field. It is easily seen that

p = p , p ⊗ p = |p ⊥ | 2 2 (I -b(x) ⊗ b(x)) + |p | 2 b(x) ⊗ b(x). Taking into account that t ∂ x b b = 1 2 ∇ x |b| 2 = 0 we deduce that p m • ∇ x f i = b(x)⊗b(x) p m •∇ x g i - (p • b(x)) |p ∧ b(x)| 2m div x b ∂ r g i + |p ∧ b(x)| 2 2m div x b ∂ z g i . (30)
and

eE(t) • ∇ p f i = e(b(x) • E(t, x)) ∂ z g i . (31)
Combining (28), ( 29), (30), (31) yields the following Vlasov equation in the phase space

(x, r, z) ∈ R 3 × R ⋆ + × R ∂ t g i + z m b(x) • ∇ x g i - zr 2m div x b ∂ r g i + r 2 2m div x b + e(b(x) • E(t, x)) ∂ z g i = 0. (32)
Notice that the magnetic momentum r 2 /(2mB(x)) is left invariant by (32). In particular (32) (supplemented by initial condition) is well posed in the phase-space (x, r, z) ∈ R 3 × R ⋆ + × R without any boundary condition at r = 0. It is possible to reformulate this equation in order to write a Vlasov equation for the dominant ion distribution f i in the phase space (x, p). For this it is sufficient to express the derivatives of g i with respect to the derivatives of f i

∂ t g i = ∂ t f i , ∂ z g i = b(x) • ∇ p f i , ∂ r g i = p -(p • b)b |p ∧ b| • ∇ p f i ∇ x g i = ∇ x f i -( ⊥ p • ∇ p f i ) t ∂ x b p |p ∧ b|
leading to the ion Vlasov equation ( 16).

Certainly the above arguments are formal. The rigorous derivation of the ion limit model ( 16) was postponed to Appendix, cf. Proposition 6.1.

Remark 3.1 As before, the invariance of the magnetic momentum |p∧b(x)| 2 /(2mB(x))

guarantees the well-posedness of (16) for p ∧ b(x) = 0. Notice also that for any (x, p)

such that p ∧ b(x) = 0, ω i (x, p) remains bounded, since ∂ x b b(x) • p |p ∧ b(x)| = ∂ x b b(x) • p -(p • b(x)) b |p ∧ b(x)| . Remark 3.2
The previous computations show that for any function f satisfying the constraint T f = 0 we have

p m • ∇ x f + eE • ∇ p f = b(x) ⊗ b(x) p m • ∇ x f + (eb ⊗ bE + ω i (x, p) ⊥ p) • ∇ p f (33)
which means that, by averaging, the transport operator p m •∇ x +eE•∇ p reduces to another transport operator, associated with the vector field (b(x)⊗b(x) p m , (eb⊗bE+ω i (x, p) ⊥ p)). An equivalent method for determining the effective transport operator is to search for a field η = (η x , η p ) such that the equality

p m • ∇ x f + eE • ∇ p f = η • ∇ x,p f (34)
holds true when f belongs to a complete family of prime integrals for T . Since in our case, (34) should be satisfied only for functions in ker T , we can assume without loss of generality that η

• (0, ω c p ∧ b) = 0, which is equivalent to η p • (p ∧ b) = 0.
Other five equations can be obtained by using the invariants x, |p ∧ b|, (p • b). Indeed, taking f = x i , i ∈ {1, 2, 3} in (34), yields

η x i = p i m = (p • b) m b i , i ∈ {1, 2, 3}. Taking f = (p • b) implies η x • ( t ∂ x b p) + η p • b = p m • ( t ∂ x b p) + eE • b = |p ∧ b| 2 2m div x b + eE • b.
Eventually one can get the last equation appealing to the invariant |p ∧ b|. Actually the computations simplify a little bit when using the invariant |p| 2 instead of |p ∧ b|. We obtain

η p • p = eE • p = e(E • b)(p • b).
Finally we retrieve the effective transport operator in (33).

Remark 3.3 It is easily seen that any prime integral of p m • ∇ x + eE • ∇ p is also a prime integral of its averaged transport operator. For example, when the electric field derives from a potential E = -∇ x φ then

p m • ∇ x + eE • ∇ p |p| 2 2m + eφ = 0 which implies b(x) ⊗ b(x) p m • ∇ x + (eb ⊗ bE + ω i (x, p) ⊥ p) • ∇ p |p| 2 2m + eφ = 0.
Remark 3.4 Nearly the same arguments apply for models with a time dependent magnetic field. In this case we have

f i (t, x, p) = g i (t, x, |p∧b(t, x)|, (p•b(t, x))
) and therefore

∂ t f i = ∂ t g i -∂ r g i (p • b) |p ∧ b| (p • ∂ t b) + ∂ z g i (p • ∂ t b).
Averaging at any fixed time t yields

∂ t f i = ∂ t g i since (p • ∂ t b) = 0.
Therefore the limit model satisfied by g i in the phase-space (x, r, z) does not change, but now, coming back in the phase-space (x, p) gives

∂ t f i + b ⊗ b p m • ∇ x f i + (eb ⊗ bE + ω i (t, x, p) ⊥ p) • ∇ p f i = 0
where, in this case, the frequency ω i is time dependent

ω i (t, x, p) = |p ∧ b(t, x)| 2m div x b - p |p ∧ b(t, x)| • ∂ t b + (∂ x b b) ⊗ b p m .

Electron limit model

In this section we derive the limit model satisfied by the dominant electron distribution f e in [START_REF] Degond | Transport coefficients of plasmas and disparate mass binary gases[END_REF] and given in Proposition 2.1 (ii). First we assume that all particles (ions and electrons) have typical momentum of the same order, that means that we are starting from the model ( 7), [START_REF] Boozer | Physics of magnetically confined plasmas[END_REF]. Next we investigate the case of comparable velocity units.

But first, one more average operator need to be introduced.

Second average operator

As has been noticed before, the analysis of the electron distribution is more complex, the Vlasov equation ( 8) involving not only the scale ε but also ε 2 . Plugging the Ansatz ( 12) into (8) leads to the time evolution equation

∂ t f e -e E • ∇ p f e + p m • ∇ x f 1 e -T f 2 e = 0 , (35)
for the dominant term f e depending on the first and second order fluctuation terms f 1 e , f 2 e , which have to be removed by using the constraints ( 13), [START_REF] Frénod | The finite Larmor radius approximation[END_REF]. The divergence constraint [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF] yields that there is a function g e = g e (t, x, r, z) such that

f e (t, x, p) = g e (t, x, |p ∧ b|, p • b).
But f e also satisfies the second constraint [START_REF] Frénod | The finite Larmor radius approximation[END_REF], given by

p m • ∇ x f e = T f 1 e = 0.
Performing the same computations as for the ion distribution f i (which is possible since the electron distribution f e also belongs to ker T ), this last constraint writes

z m b(x) • ∇ x g e - zr 2m div x b ∂ r g e + r 2 2m div x b ∂ z g e = 0. (36)
Therefore the second average operator to be considered, in order to eliminate the fluctuation term f 1 e of (35), will be that with respect to the characteristic flow

dX ds = Z(s) m b(X(s)) (37)
dR ds = - Z(s)R(s) 2m div x b (X(s)) (38) dZ ds = R(s) 2 2m div x b (X(s)) (39) (X, R, Z)(0) = (x, r, z). (40)
Observe that r 2 + z 2 is a prime integral for the field

z m b(x) • ∇ x - zr 2m div x b ∂ r + r 2 2m div x b ∂ z ( 41
)
and therefore R 2 (s) + Z 2 (s) is left invariant by the flow (37), (38), (39). The magnetic momentum

µ(x, p) = |p ∧ b(x)| 2 2mB(x) = r 2 2mB(x)
is another invariant, provided that the magnetic field is divergence free. Notice that the characteristic flow (37), (38), (39) preserves the measure dν = 2πrdxdrdz i.e.,

R 3 ×R + ×R χ((X, R, Z)(s; x, r, z)) dν = R 3 ×R + ×R χ(x, r, z) dν , ∀s ∈ R , for any function χ ∈ L 1 (R 3 × R + × R ; dν). Indeed, for any χ ∈ L 1 (R 3 × R + × R ; dν)
the solution of (32) corresponding to the initial condition χ is given by g(t, x, r, z) = χ((X, R, Z)(-t; x, r, z)).

Observe that (32) can be written into conservative form

∂ t (2πrg) + div x 2πrg z m b -∂ r 2πrg zr 2m div x b + ∂ z 2πrg r 2 2m div x b = 0 implying that d dt R 3 ×R + ×R g(t, x, r, z) dν = 0, t ∈ R
and therefore

R 3 ×R + ×R χ((X, R, Z)(-t; x, r, z)) dν = R 3 ×R + ×R χ(x, r, z) dν, t ∈ R.
We need two other invariants for solving (37), ( 38), (39). Generally a confinement region is supposed to be filled by nested magnetic surfaces, each surface enclosing the next. We consider here a simplified geometry, i.e., the framework of 2π periodic functions with respect to x 3 and cylindrical magnetic surfaces with axis parallel to e 3 = (0, 0, 1). More precisely assume that the field of unitary vectors b is given by b

(x) = ⊥ x + e 3 (x 2 1 + x 2 2 + 1) 1/2 , ⊥ x = (x 2 , -x 1 , 0). (42)
It is easily seen, by direct computations, that the scalar functions B, which are 2π periodic with respect to x 3 and satisfy the constraint div x (Bb) = 0, are those depending only on ρ and x 3 + θ, where

x 1 = ρ cos θ, x 2 = ρ sin θ. Indeed, since div x b = 0 the divergence constraint div x (Bb) = 0 is equivalent to b • ∇ x B = 0, that is x 2 ∂ x 1 B -x 1 ∂ x 2 B + ∂ x 3 B = 0.
Taking into account that

∂ θ = x 1 ∂ x 2 -x 2 ∂ x 1 ,
we obtain

-∂ θ B + ∂ x 3 B = 0 ,
whose invariants are ρ = (x 2 1 + x 2 2 ) 1/2 and x 3 + θ.

Proposition 4.1 Assume that the direction of the magnetic field is given by (42). The characteristic flow of (37), (38), (39) is then given by

X 1 (s; x, r, z) = x 1 cos sz m 1 + ρ 2 + x 2 sin sz m 1 + ρ 2 , X 2 (s; x, r, z) = -x 1 sin sz m 1 + ρ 2 + x 2 cos sz m 1 + ρ 2 , X 3 (s; x, r, z) = x 3 + sz m 1 + ρ 2 , R(s; x, r, z) = r, Z(s; x, r, z) = z. Proof Since div x b = 0, we have (R, Z)(s; x, r, z) = (r, z). The components X 1 , X 2 , X 3 satisfies dX 1 ds = z m X 2 (s) 1 + X 1 (s) 2 + X 2 (s) 2 , dX 2 ds = - z m X 1 (s) 1 + X 1 (s) 2 + X 2 (s) 2 , dX 3 ds = z m 1 1 + X 1 (s) 2 + X 2 (s) 2 .
Clearly ρ = x 2 1 + x 2 2 is left invariant and therefore our conclusion follows immediately. Notice that for any z = 0 the characteristics (X, R, Z)(s; x, r, z) are T (ρ, z) = 2π z m 1 + ρ 2 periodic in R 2 × R/(2πZ) × R + × R and that for z = 0 we have (X, R, Z)(s; x, r, 0) = (x, r, 0).

Let us now introduce the average operator with respect to the flow given in Proposition 4.1. It is easily seen that the 2π periodic functions with respect to x 3 and constant along the above flow are those depending only on ρ, x 3 + θ, r and z. We introduce the first order differential operator

T 1 v = div x z m ⊥ x + e 3 1 + x 2 1 + x 2 2 v , D(T 1 ) = v ∈ L 2 # (R 3 × R + × R ; dν) : div x z m ⊥ x + e 3 1 + x 2 1 + x 2 2 v ∈ L 2 # (R 3 × R + × R ; dν) .
Here L 2 # (R 3 × R + × R ; dν) stands for the space of 2π periodic functions with respect to x 3 , measurable and such that

v 1 := R 2 2π 0 R + R |v(x, r, z)| 2 dν 1/2
< +∞.

To this differential operator we associate the average operator along the trajectories of Proposition 4.1 (see [START_REF] Bostan | Transport equations with singular coefficients. Application to the gyrokinetic models in plasma physics[END_REF]), given by

v 1 (x, r, z) = lim T →+∞ 1 T T 0 v((X, R, Z)(s; x, r, z)) ds.
If z = 0 (and thus for a negligible set) one gets v 1 (x, r, 0) = v(x, r, 0) and if z = 0 we have

v 1 (x, r, z) = 1 T (ρ, z) T (ρ,z) 0 v((X, R, Z)(s; x, r, z)) ds = 1 2π 2π 0 v(R(-α) t (x 1 , x 2 ), x 3 + α, r, z) dα = 1 2π 2π 0 v(R(α) t (x 1 , x 2 ), x 3 -α, r, z) dα where R(α) is the rotation matrix of angle α R(α) =   cos α -sin α sin α cos α   , α ∈ R.
Actually v 1 depends only on ρ, x 3 + θ, r, z

Proposition 4.2

The average operator • 1 is linear and continuous on L 2 # (R 3 × R + × R ; dν). Moreover it coincides with the orthogonal projection on the kernel of T 1 i.e.,

• 1 : L 2 # (R 3 × R + × R ; dν) → ker T 1 and R 2 2π 0 R + R (v -v 1 )ψ dν = 0, ∀ψ ∈ ker T 1 . Proof For any function v ∈ L 2 # (R 3 × R + × R ; dν) we have | v 1 | 2 (x, r, z) ≤ 1 2π 2π 0 v 2 (ρ cos α, ρ sin α, x 3 + θ -α, r, z) dα.
Integrating with respect to x ∈ R 2 × [0, 2π[ and using polar coordinates one gets for almost all (r, z)

∈ R + × R R 2 2π 0 | v 1 | 2 dx ≤ 1 2π 2π 0 R + 2π 0 2π 0 v 2 (ρ cos α, ρ sin α, x 3 + θ -α, r, z)ρ dx 3 dθdρ dα = R 2 2π 0 v 2 (x, r, z) dx.
Multiplying by 2πr and integrating with respect to (r,

z) ∈ R + × R yields R 2 2π 0 R + R v 2 1 dν ≤ R 2 2π 0 R + R v 2 2πr dzdrdx = R 2 2π 0 R + R v 2 dν saying that • 1 ∈ L(L 2 # (R 3 × R + × R ; dν), L 2 # (R 3 × R + × R ; dν)) and • 1 L ≤ 1. Pick a function ψ ∈ ker T 1 i.e., ∃ χ : ψ(x, r, z) = χ(ρ, x 3 + θ, r, z) ∈ L 2 # (R 3 × R + × R ; dν) and let us show that R 2 2π 0 R + R (v -v 1 )ψ dν = 0.
We are done if we prove that Taking the average with respect to α ∈ [0, 2π[ we obtain

R 2 2π 0 (v -v 1 )ψ dx = 0, (r, z) ∈ R + × R. For (r, z, α) ∈ R + × R × [0,
R 2 2π 0 vψ dx = R + 2π 0 2π 0 χ(ρ, x 3 + θ, r, z) 1 2π 2π 0 v(ρ cos(θ -α), ρ sin(θ -α), x 3 + α, r, z) dαρdx 3 dθdρ = R + 2π 0 2π 0 χ(ρ, x 3 + θ, r, z) v 1 (ρ cos θ, ρ sin θ, x 3 , r, z)ρ dx 3 dθdρ = R 2 2π 0 v 1 ψ dx since 1 2π 2π 0 v(ρ cos(θ -α), ρ sin(θ -α), x 3 + α, r, z) dα = 1 2π θ θ-2π v(ρ cos u, ρ sin u, x 3 + θ -u, r, z) du = 1 2π 2π 0 v(ρ cos u, ρ sin u, x 3 + θ -u, r, z) du = v 1 (ρ cos θ, ρ sin θ, x 3 , r, z).

The previous computations show that

R 2 2π 0 R + R (v-v 1 )ψ dν = 0 for any ψ ∈ ker T 1 saying that v 1 = Proj ker T 1 v for any v ∈ L 2 # (R 3 × R + × R ; dν). In particular v 1 = v for any v ∈ ker T 1 and • 1 L = 1 (notice that ker T 1 = ∅, for example exp(-x 2 1 -x 2 2 -r 2 -z 2 ) ∈ ker T 1 ). Remark 4.1
The key point in the construction of the average operator • 1 is that the measure ν is left invariant by the flow (37), ( 38), (39). The reader can convince himself that for any function ψ ∈ ker T 1 the following formal computations hold true (see also [START_REF] Bostan | Transport equations with singular coefficients. Application to the gyrokinetic models in plasma physics[END_REF], [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF] pp. 57)

R 2 2π 0 R + R v 1 ψ dν = lim T →+∞ 1 T T 0 R 2 2π 0 R + R v((X, R, Z)(s; x, r, z))ψ(x, r, z) dν dt = lim T →+∞ 1 T T 0 R 2 2π 0 R + R (v ψ)((X, R, Z)(s; x, r, z)) dν dt = lim T →+∞ 1 T T 0 R 2 2π 0 R + R v(x, r, z)ψ(x, r, z) dν dt = R 2 2π 0 R + R v(x, r, z)ψ(x, r, z) dν.
The previous result also gives the orthogonal decomposition of L 2 # (R 3 × R + × R ; dν) into invariant functions along the flow (37), (38), (39) and zero average functions i.e.,

L 2 # (R 3 × R + × R ; dν) = ker T 1 ⊕ ⊥ ker • 1 since v = v 1 + (v -v 1 ), R 2 2π 0 R + R (v -v 1 ) v 1 dν = 0. Observing that T 1 is skew-adjoint on L 2 # (R 3 × R + × R ; dν) we also deduce that ker • 1 = ker(Proj ker T 1 ) = (ker T 1 ) ⊥ = (ker T ⋆ 1 ) ⊥ = Range T 1 .
In particular Range T 1 ⊂ ker • 1 .

Comparable ion/electron momentum units

The derivation of the electron limit model requires long computations, since double average is needed. For the sake of the presentation, it is done in several steps.

Proof of Proposition 2.1 (ii)

Let us assume that the magnetic field is given by (42). Recall that f e satisfies the constraint [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF] saying that

f e (t, x, p) = g e (t, x, r = |p ∧ b(x)|, z = p • b(x) ) ,
and that g e verifies the constraint (36) (which is a consequence of ( 14)) implying that

g e (t, x, r, z) = h e (t, ρ = (x 2 1 + x 2 2 ) 1/2 , s = x 3 + θ, r, z).
Elimination of the second order correction f 2 e Before investigating (35) observe that ( 14) allows us to determine the zero average part of the first order correction f 1 e in terms of the leading order distribution f e . Indeed we have the decomposition in

L 2 (R 3 × R 3 ) f 1 e = g 1 e (t, x, r = |p ∧ b|, z = p • b ) + f 1 e , f 1 e = 0
and therefore by ( 14) one gets

T f 1 e = p m • ∇ x f e , f 1 e = 0. ( 43
)
By (36) we know that p m • ∇ x f e = 0 and thus Proposition 3.2 guarantees the solvability of (43). Actually we can write

p m • ∇ x f e = p m • ∇ x f e - p m • ∇ x f e = p m - p m • ∇ x g e -∂ r g e p • b |p ∧ b| ∂ x b : p ⊗ p -p ⊗ p m + ∂ z g e ∂ x b : p ⊗ p -p ⊗ p m
and straightforward computations imply

T -1 (p -p ) = - p ∧ b ω c (44) T -1 (p ⊗ p -p ⊗ p ) = - p ∧ b ω c ⊗ 3 4 (p • b) b + p 4 - 3 4 (p • b) b + p 4 ⊗ p ∧ b ω c . ( 45
)
Therefore the zero average distribution f 1 e is given by

f 1 e = - p ∧ b mω c • ∇ x g e -∂ z g e -∂ r g e p • b |p ∧ b| ∂ x b + t ∂ x b : p ∧ b mω c ⊗ 3 4 (p • b) b + p 4 . ( 46
)
The time evolution equation for f e comes by (35) after eliminating the distributions f 1 e , f 2 e . Applying the average operator • allows us to get rid of

f 2 e ∂ t f e -eE • ∇ p f e + p m • ∇ x g 1 e (t, x, |p ∧ b|, p • b) + p m • ∇ x f 1 e = 0. ( 47
)
Since f e (t, x, p) and g 1 e (t, x, |p ∧ b|, p • b ) satisfy the constraint [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF] we have as in ( 29), (30), (31) (and by taking into account that div x b = 0)

∂ t f e = ∂ t g e , e E • ∇ p f e = e(b • E) ∂ z g e , p m • ∇ x g 1 e = b ⊗ b p m • ∇ x g 1 e .
Plugging the above expressions into (47) yields

∂ t g e -e (b • E) ∂ z g e + p m • ∇ x f 1 e + T 1 g 1 e = 0. ( 48
)
Applying • 1 to this equation will eliminate the term T 1 g 1 e . The difficult task is now to give an explicit formula for the average p m • ∇ x f 1 e in terms of g e by using (43). It happens that this average also reduces to a transport operator ξ • ∇ x,r,z in the phase space (x, r, z) (see Remark 3.2).

Computation of the field ξ

Indeed, for any function g e = g e (x, r, z) satisfying the constraint T 1 g e = 0 we have

p m • ∇ x f 1 e = ξ • ∇ x,r,z g e (49)
where f 1 e is the unique solution of

f 1 e = 0, T f 1 e = p m • ∇ x g e + ∂ z g e - p • b |p ∧ b| ∂ r g e ∂ x b : p ⊗ p m .
Since (49) has to be satisfied only for functions g e ∈ ker T 1 we can assume that the field ξ = (ξ x , ξ r , ξ z ) verifies

ξ x • b = 0. ( 50
)
Other four equalities are obtained by imposing (49) when g e is one of the invariants

ρ = x 2 1 + x 2 2 , s = x 3 + θ, r, z.
For example taking g e = (r 2 + z 2 )/2 we get f 1 e = 0 and thus rξ r + zξ z = 0.

(51)

Consider now g e = ρ. In this case (46) gives

f 1 e = - ∇ x ρ m • p ∧ b ω c = p mω c • (∇ x ρ ∧ b)
and therefore

p m • ∇ x f 1 e = - p m 2 ω c • ∇ x B B p • (∇ x ρ ∧ b) + p m 2 ω c • t ∂ x (∇ x ρ ∧ b)p. ( 52
)
By direct computations one gets

∇ x ρ ∧ b = 1 ρ b - 1 + ρ 2 ρ e 3 ∂ x (∇ x ρ ∧ b) = 1 ρ ∂ x b - 1 ρ 2 b ⊗ ∇ x ρ -e 3 ⊗ ∇ x 1 + ρ 2 ρ Tr(∂ x (∇ x ρ ∧ b)) = 1 ρ div x b - 1 ρ 2 b • ∇ x ρ -∂ x 3 1 + ρ 2 ρ = 0 b • (∂ x (∇ x ρ ∧ b) b) = 1 ρ b • ∂ x b b - 1 ρ 2 ∇ x ρ • b -b 3 b • ∇ x 1 + ρ 2 ρ = 0.
Notice that for any matrix A(x) and vectors η(x), χ(x) we have

p ⊗ p : A(x) = |p ∧ b| 2 2 (Tr(A) -b • Ab) + (p • b) 2 b • Ab (53) (p • η)(p • χ) = |p ∧ b| 2 2 (η • χ -(η • b)(χ • b)) + (p • b) 2 (η • b)(χ • b). ( 54
)
Taking the average of (52) yields

p m • ∇ x f 1 e = - |p ∧ b| 2 2m 2 ω 2 c ∇ x ω c • (∇ x ρ ∧ b) = |p ∧ b| 2 2m 2 ω c 1 + ρ 2 ρ ∂ s B B
and therefore (49) implies

ξ x 1 x 1 ρ + ξ x 2 x 2 ρ = r 2 2m 2 ω c 1 + ρ 2 ρ ∂ s B B . (55)
Consider now g e = s = x 3 + θ. In this case

f 1 e = 0, T f 1 e = p 3 m + p 2 x 1 -p 1 x 2 mρ 2
and therefore

f 1 e = - (p ∧ b) 3 mω c + ⊥ x mω c ρ 2 • (p ∧ b) = p 1 x 1 + p 2 x 2 mω c ρ 2 1 + ρ 2 .
We deduce that

p m •∇ x f 1 e = p m • (p 1 , p 2 , 0) mω c ρ 2 1 + ρ 2 + p 1 x 1 + p 2 x 2 mω c ∇ x 1 + ρ 2 ρ 2 - ∇ x B B 1 + ρ 2 ρ 2 implying that p m • ∇ x f 1 e = 1 + ρ 2 m 2 ω c ρ 2 p 2 1 + p 2 2 - |p ∧ b| 2 2m 2 ω c 2 + ρ 2 ρ 2 1 + ρ 2 + 1 + ρ 2 ρ 2 x 1 ∂ x 1 B + x 2 ∂ x 2 B B = (p • b) 2 m 2 ω c 1 + ρ 2 - |p ∧ b| 2 2m 2 ω c ∂ ρ B B 1 + ρ 2 ρ .
Consequently the choice g e = s = x 3 + θ in (49) leads to

- x 2 ρ 2 ξ x 1 + x 1 ρ 2 ξ x 2 + ξ x 3 = 1 m 2 ω c 1 + ρ 2 z 2 - ∂ ρ B B 1 + ρ 2 ρ r 2 2 . ( 56
)
Finally taking g e = z we obtain

T f 1 e = ∂ x b : p ⊗ p m , f 1 e = 0
and by ( 45)

f 1 e = 1 4mω c (I + 3b ⊗ b)(∂ x b + t ∂ x b )M [b] : p ⊗ p where M [b] is the matrix of the linear map p → b ∧ p M [b] =      0 -b 3 b 2 b 3 0 -b 1 -b 2 b 1 0     
.

By direct computations one gets

1 4 (I + 3b ⊗ b)(∂ x b + t ∂ x b )M [b] = b ⊗ e 3 1 + ρ 2 - b ⊗ b 1 + ρ 2
and thus the distribution f 1 e can be written

f 1 e = 1 mω c p 3 (p • b) 1 + ρ 2 - (p • b) 2 1 + ρ 2 . The transport term p m • ∇ x f 1 e appears like p m • ∇ x f 1 e = - p m 2 ω c • ∇ x B B p 3 (p • b) 1 + ρ 2 - (p • b) 2 1 + ρ 2 + p m 2 ω c d dρ 1 1 + ρ 2 (p • b) b 3 ∇ x ρ - d dρ 1 1 + ρ 2 (p • b) 2 ∇ x ρ + p m 2 ω c • p 3 1 + ρ 2 - 2(p • b) 1 + ρ 2 t ∂ x b p. (57)
It is easily seen that

p • ∇ x B p 3 = |p ∧ b| 2 2 ∂ s B. p • ∇ x B = (p • b)(b • ∇ x B) = 0, p • ∇ x ρ = (p • b)(b • ∇ x ρ) = 0 p • t ∂ x b p = ∂ x b : p ⊗ p = 0.
For the last term in (57) use the formula for k ∈ {1, 2, 3}

p ⊗ p p k = (p • b)b k (p • b) 2 - 3 2 |p ∧ b| 2 b ⊗ b + (p • b)|p ∧ b| 2 2 [b k I + b ⊗ e k + e k ⊗ b]
to obtain

p 3 (p • t ∂ x b p) = (p • b)|p ∧ b| 2 2 (∂ x b b) 3 + ( t ∂ x b b) 3 = 0.
Taking the average of (57) we deduce

p m • ∇ x f 1 e = - (p • b)|p ∧ b| 2 2m 2 ω c 1 + ρ 2 ∂ s B B
and therefore the equality (49) with g e = z becomes

ξ z = - r 2 z 2m 2 ω c 1 + ρ 2 ∂ s B B . (58)
The solution of the linear system (50), (51), (55), (56), (58) is

ξ x = z 2 m 2 ω c (∂ x b b ∧ b) - r 2 2m 2 ω c b ∧ ∇ x B B ξ r = r z 2 2m 2 ω c ∇ x B B • (∂ x b b ∧ b), ξ z = - r 2 z 2m 2 ω c ∇ x B B • (∂ x b b ∧ b).
Remark 4.2 Notice that it is also possible to compute the field ξ in (49) by working on p m • ∇ x f 1 e , with f 1 e coming from (46), but the computations would be much more complex. Actually we have used the invariants of T 1 only for simplifying the calculations. Therefore, up to this step we don't need a complete family of invariants to be available for T 1 . The computations can be done for general magnetic shapes (recall that the main motivation when restricting to [START_REF] Hazeltine | Plasma confinement[END_REF] was the need of invariants).

Elimination of the first order correction

∂ t g e -e (b • E) ∂ z g e + ξ • ∇ x,r,z g e + T 1 g 1 e = 0 (59)
under the constraint T 1 g e = 0. The last step to be accomplished is to eliminate T 1 g 1 e in (59) by applying the average operator • 1

∂ t g e -e (b • E) ∂ z g e + ξ • ∇ x,r,z g e 1 = 0, T 1 g e = 0.
We obtain thus

∂ t g e + χ • ∇ x,r,z g e = 0 (60)
for some field χ = (χ x , χ r , χ z ) satisfying

-e (b • E) ∂ z g e + ξ • ∇ x,r,z g e 1 = χ • ∇ x,r,z g e (61)
for any function g e ∈ ker T 1 (see Remark 3.2).

Computation of the field χ

It is sufficient to impose (61) when g e belongs to the family of invariants ρ =

x 2 1 + x 2 2 , s = x 3 + θ, r, z and to assume that χ x • b = 0. Taking g e = r and g e = z one gets

χ r = ξ r 1 , χ z = -e (b • E) + ξ z 1 .
Since B depends only on the invariants ρ, s it is easily seen that ξ r 1 = ξ r ,

ξ z 1 = ξ z and b • E(t) 1 = 1 2π 2π 0 (b • E(t))(ρ cos α, ρ sin α, s -α) dα.
Assuming that E(t, x) = -∇ x φ(t, x) for some 2π periodic potential with respect to x 3 and observing that

d dα {φ(t, ρ cos α, ρ sin α, s -α)} = 1 + ρ 2 (E • b)(ρ cos α, ρ sin α, s -α) we deduce that b • E(t) 1 = 0.
Applying now (61) with g e = ρ leads to

χ x 1 x 1 ρ + χ x 2 x 2 ρ = ξ x 1 x 1 ρ + ξ x 2 x 2 ρ 1 = r 2 2m 2 ω c 1 + ρ 2 ρ ∂ s B B .
Finally the choice g e = s = x 3 + θ implies

-χ x 1 x 2 ρ 2 + χ x 2 x 1 ρ 2 + χ x 3 = -ξ x 1 x 2 ρ 2 + ξ x 2 x 1 ρ 2 + ξ x 3 1 = 1 m 2 ω c 1 + ρ 2 z 2 - r 2 (1 + ρ 2 ) 2ρ ∂ ρ B B .
We have obtained the same equations as for ξ (this is due to the fact that the magnetic intensity depends only on ρ, s and that the electric field derives from a potential, which

ensures that b • E(t) 1 = 0). Therefore χ = ξ.
Remark 4.3 For computing the field χ in (61) we really need a complete family of invariants for T 1 . Therefore we have to restrict ourselves to particular magnetic shapes, for example [START_REF] Hazeltine | Plasma confinement[END_REF].

Vlasov equation for the leading order electron density f e

Let us now re-write the equation (60) in the standard phase space (x, p) by using the formulae

∂ t g e = ∂ t f e , ∂ z g e = b•∇ p f e , ∂ r g e = p -(p • b) b |p ∧ b| •∇ p f e , ∇ x g e = ∇ x f e -⊥ p•∇ p f e t ∂ x b p |p ∧ b| .
Thus we obtain the electron Vlasov equation

∂ t f e + (v GD + v CD ) • ∇ x f e + ω e (x, p) ⊥ p • ∇ p f e = 0 (62)
where v GD (resp. v CD ) are the magnetic gradient (resp. curvature ) electron drifts

v GD = |p ∧ b| 2 2m 2 (-ω c ) b ∧ ∇ x B B , v CD = (p • b) 2 m 2 (-ω c ) (b ∧ ∂ x b b)
and the frequency ω e is given by

ω e (x, p) = -(v GD + v CD ) • t ∂ x b p |p ∧ b| + |p ∧ b|(p • b) 2m 2 ω c ∇ x B B • (b ∧ ∂ x b b).
The just presented asymptotic analysis has been performed under the assumption that all particles (ions and electrons) have typical momentum of the same order, saying that the electron velocity is much larger than the ion velocity

v e v i = m i m e = 1 µ >> 1.
Averaging with respect to the fast cyclotronic motion led to the ion model

∂ t f i + b(x) ⊗ b(x) p m • ∇ x f i + eb(x) ⊗ b(x)E + ω i (x, p) ⊥ p • ∇ p f i = 0 , (63)
whereas the double averaging yield the Vlasov equation (62) for the electron motion.

The behaviour of ions and electrons are very different; the ions are advected along the parallel direction whereas the electrons are advected along the orthogonal directions with respect to the magnetic field. Only the ions remain confined at the leading order around the magnetic lines, whereas the electrons are submitted to orthogonal drifts (magnetic gradient drift, magnetic curvature drift) at the leading order and not at the next order, as for the ions. This is precisely due to the assumption that the typical electron velocity is much larger than the typical ion velocity.

Comparable ion/electron velocity units

We will investigate the asymptotic behaviour of ( 21), [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF] for general smooth magnetic shapes and constant electric potentials along the magnetic lines

E(t, x) = -∇ x φ(t, x), b • ∇ x φ = 0. ( 64
)
We will see that neglecting the parallel component (with respect to the magnetic field)

of the electric field will simplify a lot the computations. In particular we don't need to use double average when determining the electron limit model, only simple average is sufficient. From the physical point of view neglecting the parallel electric field comes by considering the MHD closure, see [START_REF] Hazeltine | Plasma confinement[END_REF] pp. 231. The fundamental assumption of the MHD consists in assuming that the velocity is given by

V = (V • b) b + E ∧ B ε |B ε | 2 . ( 65
)
Alternatively we write the MHD Ohm's law E + V ∧ B ε = 0 thus ensuring both (65)

and

E • b = 0.
Obviously the study of ( 21) is identical to that of (7) leading to a limit model similar to [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF] but in the phase space (x, v)

∂ t F i + b(x) ⊗ b(x)v • ∇ x F i + Ω i (x, p) ⊥ v • ∇ v F i = 0 ( 66
)
where

Ω i (x, p) = |v ∧ b(x)| 2 div x b(x) -(v • b(x)) ∂ x b b(x) • v |v ∧ b(x)| and ⊥ v = |v ∧ b(x)| b(x) -(v • b(x)) v -(v • b(x)) b(x) |v ∧ b(x)| .
It remains to analyze [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF]. Plugging the Ansatz F ε e = F e + εF 1 e + ε 2 F 2 e + ... in ( 22) yields

T F e = ω c (x)(v ∧ b) • ∇ v F e = 0 ( 67
) - eE m • ∇ v F e -T F 1 e = 0 ( 68
)
∂ t F e + v • ∇ x F e - eE m • ∇ v F 1 e -T F 2 e = 0 (69)
. . . The constraint T F e = 0 implies that there is a function G e = G e (t, x, r, z) depending on time t and the invariants x, r

= |v ∧ b(x)|, z = v • b(x) such that F e (t, x, v) = G e (t, x, |v ∧ b(x)|, v • b(x) ).
As before let us decompose the density F 1 e into a constant part along the flow of T and a zero average part

F 1 e (t, x, v) = G 1 e (t, x, r = |v ∧ b(x)|, z = v • b(x) ) + F 1 e , F 1 e = 0.
Observe that by (68), (64) we can write

T F 1 e = - eE m • v -(v • b)b |v ∧ b| ∂ r G e = eE mω c |v ∧ b| • T (v ∧ b) ∂ r G e (70)
which allows us to determine the zero average part of

F 1 e F 1 e = eE mω c |v ∧ b| • (v ∧ b) ∂ r G e = - (v ED • v) |v ∧ b| ∂ r G e ( 71
)
where v ED = E∧b B is the electric cross field drift. Actually the computation (70) shows that eE m • ∇ v F ∈ Range T = ker • for any function F ∈ ker T . Therefore applying the average operator to (68) doesn't yield any other constraint for the leading order term F e . In particular we have

eE m • ∇ v {G 1 e (t, x, |v ∧ b(x)|, v • b(x) )} = 0
and we don't need double averaging for identifying the time evolution equation for F e .

Indeed, taking the average of (69) allows to eliminate both eE m

• ∇ v G 1 e and T F 2 e ∂ t F e + v • ∇ x F e - eE m • ∇ v F 1 e = 0. (72)
We are done if we compute the average of eE m • ∇ v F 1 e in terms of the density F e by using (71).

Proposition 4.3 Assume that F 1 = -∂ r G(x, |v ∧ b|, v • b)(v ED • v)/|v ∧ b| for some function G = G(x, r, z). Then we have eE m • ∇ v F 1 = 0.
Proof We check easily that

∇ v {∂ r G(x, |v ∧ b|, v • b)} = ∂ 2 r G v -(v • b)b |v ∧ b| + ∂ 2 rz G b ∇ v (v ED • v) |v ∧ b| = v ED |v ∧ b| - (v ED • v) |v ∧ b| 2 v -(v • b) b |v ∧ b| and therefore -∇ v F 1 = (v ED • v) |v ∧ b| ∂ 2 r G v -(v • b)b |v ∧ b| + ∂ 2 rz G b + ∂ r G |v ∧ b| v ED - (v ED • v) |v ∧ b| v -(v • b) b |v ∧ b| . Using the formula v ⊗ v = |v ⊥ | 2 2 (I -b ⊗ b) + |v | 2 b ⊗ b yields -∇ v F 1 = ∂ 2 r G 2 v ED + ∂ r G 2|v ∧ b| v ED and finally eE m • ∇ v F 1 = 0.
Performing similar computations as those in (30) we deduce that

v • ∇ x F e = (b ⊗ b)v • ∇ x G e - (v • b)|v ∧ b| 2 div x b ∂ r G e + |v ∧ b| 2 2 div x b ∂ z G e
∂ t F e + b(x) ⊗ b(x)v • ∇ x F e + Ω i (x, p) ⊥ v • ∇ v F e = 0. (73)
Notice that in this case the leading order distributions F i /F e for the ions/electrons satisfy the same model (cf. 0) in Proposition 2.2). This is due to the fact that the electric potential is constant along the magnetic lines i.e., the electric field does not accelerate any particle along the magnetic lines. Comparing [START_REF] Rax | Physique des plasmas, Cours et applications[END_REF], [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF] we may expect a different behaviour between ions and electrons, since the electron Coulomb acceleration is much stronger than the ion Coulomb acceleration (because m e << m i ) but this occurs only in the orthogonal directions with respect to the magnetic lines. And since at the leading order (i.e., after averaging with respect to the fast cyclotronic motion) the orthogonal electric field doesn't play any role to the particle dynamics, we obtain similar Vlasov equations for ions and electrons cf. (66), (73).

Drift velocities

At the leading order, both particle species are confined along the magnetic lines. But specific drift velocities in the orthogonal directions are expected at the next order.

Indeed, let us compute the current densities of the first order corrections.

Proof of Proposition 2.2 (i) Multiplying (68) by v yields after integration

R 3 F e eE m dv + ω c R 3 F 1 e v ∧ b dv = 0 implying that R 3 F 1 e (v -(v • b) b) dv = R 3 F e dv E ∧ b B .
Taking into account that the dominant electron density has no current in the orthogonal direction we deduce that the electron mean velocity in the orthogonal directions is given by the electric cross field drift v ED as

R 3 (F e + εF 1 e )(v -(v • b) b) dv = ε R 3 F 1 e (v -(v • b) b) dv = ε R 3 F e dv v ED . ( 74
)
It remains to compute the ion drifts. Multiplying (11) (written in (x, v) phase space) by v and integrating with respect to the velocity yields

∂ t R 3 F i v dv + div x R 3 F i v ⊗ v dv - eE m R 3 F i dv -ω c R 3 F 1 i (v ∧ b) dv = 0. ( 75
)
As mentioned before the current of the dominant density F i is parallel to the magnetic field

R 3 F i v dv = R 3 F i (v • b) dv b
and thus, by taking the vector product of (75) by b, we obtain

R 3 F 1 i (v -(v • b) b) dv = R 3 F i dv E ∧ b B + b ω c ∧ div x R 3 F i v ⊗ v dv. (76)
The last term in the above equality can be expressed in terms of magnetic gradient drift v GD and magnetic curvature drift v CD .

Lemma 4.1 Let us consider the orthogonal (resp. parallel) velocity V ⊥ (resp. V )

given by Proof We have

V ⊥ (t, x) = R 3 |v ∧ b(x)| 2 F i (t, x, v) dv R 3 F i (t, x, v) dv 1/2 ∈ R + V (t, x) = R 3 (v • b(x)) 2 F i (t, x, v) dv R 3 F i (t, x, v) dv 1/2 ∈ R + . Then we have b ω c ∧ div x R 3 F i v ⊗ v dv = b ∧ ∇ x R 3 µ(x, v) e F i dv + V 2 ⊥ 2ω c b ∧ ∇ x B B R 3 F i dv + V 2 ω c - V 2 ⊥ 2ω c k(x)b(x) ∧ n(x) R 3 F i dv (77
R 3 F i v ⊗ v dv = R 3 F i (v • b) 2 dv b ⊗ b + R 3 F i |v ∧ b| 2 2 dv (I -b ⊗ b) = V 2 R 3 F i dv b(x) ⊗ b(x) + V 2 ⊥ 2 R 3 F i dv (I -b ⊗ b).
Using the formula div x (aA) = A(x)∇ x a + a(x)div x A for any smooth scalar function a(x) and matrix function A(x) and taking into account that div

x b ⊗ b = ∂ x bb + div x b b one gets b ω c ∧ div x V 2 R 3 F i dv b ⊗ b = V 2 R 3 F i dv ω c k(x)b(x) ∧ n(x).
Similarly we obtain

b ω c ∧ div x V 2 ⊥ 2 R 3 F i dv (I -b ⊗ b) = b ∧ ∇ x V 2 ⊥ 2ω c R 3 F i dv + V 2 ⊥ 2ω c R 3 F i dv b ∧ ∇ x B B - V 2 ⊥ 2ω c R 3 F i dv k(x)b(x) ∧ n(x)
and our conclusion follows immediately.

Based on Lemma 4.1 we finish the proof of Proposition 2.2.

Proof of Proposition 2.2 (ii) Combining ( 76), (77) yields

R 3 (F i + εF 1 i )(v -(v • b) b) dv = ε R 3 F 1 i (v -(v • b) b) dv = ε R 3 F i dv E ∧ b B + V 2 ⊥ 2ω c b ∧ ∇ x B B + V 2 ω c - V 2 ⊥ 2ω c k(x)b(x) ∧ n(x) + εb ∧ ∇ x R 3 µ(x, v) e F i dv
and therefore the ion mean velocity in the orthogonal directions is given by

ε E ∧ b B + V 2 ⊥ 2ω c b ∧ ∇ x B B + V 2 ω c - V 2 ⊥ 2ω c k(x)b(x) ∧ n(x) + b ∧ ∇ x R 3 µ(x, v)F i dv e R 3 F i dv . (78)
The formula (78), (74) clearly emphasize the specific behaviours of the ions and electrons. When first order corrections are taken into account the ions deviate along the orthogonal directions and both electric and magnetic drifts are observed. The electrons deviate as well along the orthogonal directions, but only the electric cross field drift occurs.

Conclusion

The aim of this paper was to rigorously investigate the ε → 0 asymptotics of the ion/electron Vlasov equations, where the small parameter ε accounts for the electron/ion mass ratio as well as the fast cyclotronic motion. Depending on the initial assumption one makes, concerning the order of magnitude of the ion/electron momenta or velocities, different limit models are obtained. Electrons and ions behave differently with regard to the considered ε-order as well as the considered motion direction (parallel or perpendicular to the magnetic field lines). The rigorous asymptotic analysis performed in this paper can be extended to several other transport problems, involving multiple scales.

Appendix

Motivated by the formal considerations in Proposition 2.1 (i) we establish here the following weak convergence result Proposition 6.1 Assume that the electro-magnetic field (E, B b) is smooth i.e., E ∈

L 1 loc (R + ; W 1,∞ (R 3 )) 3 , b ∈ W 2,∞ (R 3 ) 3 , B ∈ W 1,∞ (R 3 ) such that inf x∈R 3 B(x) > 0. For any ε > 0 let f ε i ∈ L ∞ (R + ; L 2 (R 3 × R 3
)) be the weak solution (by characteristics) of (7) satisfying the initial condition

f ε i (0) = f in i ∈ L 2 (R 3 × R 3 ). Then there is a sequence ε n ց 0 such that (f εn i ) n converges weakly ⋆ in L ∞ (R + ; L 2 (R 3 × R 3
)) towards the weak solution of the problem [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF], understood in D ′ ({(t, x, p) : p ∧ b(x) = 0}), with the initial condition f i (0) = f in i .

Proof Clearly we have for any ε > 0, t ∈ R +

R 3 R 3 (f ε i (t, x, p)) 2 dpdx = R 3 R 3 (f in i (x, p)) 2 dpdx
and therefore there is a sequence ε n ց 0 and f i ∈ L ∞ (R + ; L 2 (R 3 × R 3 )) such that lim n→+∞ f εn i = f i weakly ⋆ in L ∞ (R + ; L 2 (R 3 × R 3 )). Using the weak formulation of [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] with test functions η(t)ϕ(x, p), η ∈ C 1 c (R + ), ϕ ∈ C 1 c (R 3 × R 3 ), we deduce after multiplication by ε n

R + η R 3 R 3 f i ω c (p ∧ b) • ∇ p ϕ dpdxdt = lim n→+∞ R + η R 3 R 3
R 3 R 3 f in i ϕ dpdx - R + η ′ (t) R 3 R 3 f εn i ϕ dpdxdt - R + η(t) R 3 R 3 f εn i p m • ∇ x ϕ + eE • ∇ p ϕ dpdxdt = 0
and letting n → +∞ yields -η(0)

R 3 R 3 f in i ϕ dpdx - R + η ′ (t) R 3 R 3 f i ϕ dpdxdt (79) - R + η(t) R 3 R 3 f i p m • ∇ x ϕ + eE • ∇ p ϕ dpdxdt = 0.
Let us denote by g in i = g in i (x, r, z) the function such that f in i (x, p) = g in i (x, |p ∧ b(x)|, (p • b(x))). Since T ϕ = 0 we have with the notation dν = 2πrdxdrdz

R 3 R 3 f in i ϕ dpdx = R 3 R 3 f in i ϕ dpdx = R 3 ×R + ×R g in i ψ dν. ( 80
)
Notice also that

R 3 R 3 f i (t,
Similarly, since T f i = 0 we have

R 3 R 3 f i p m • ∇ x ϕ + eE • ∇ p ϕ dpdx = R 3 R 3 f i p m • ∇ x ϕ + eE • ∇ p ϕ dpdx.
Notice that ϕ is smooth, vanishes in a neighborhood of the set {(x, p) : p ∧ b(x) = 0}

and T ϕ = 0. Therefore we can apply (30), (31), leading to the formula

R 3 R 3 f i p m • ∇ x ϕ + eE • ∇ p ϕ dpdx (82) = R 3 ×R + ×R g i z m b • ∇ x ψ - zr 2m div x b ∂ r ψ + r 2 2m div x b + e(b • E) ∂ z ψ dν.
Combining (79), (80), (81), (82) yields η(0) -η(0)

R 3 ×R + ×R g in i ψ dν - R + η ′ (t) R 3 ×R + ×R g i ψ dνdt (83) - R + η R 3 ×R + ×R g i z m b • ∇ x ψ - zr 2m div x b ∂ r ψ +
R 3 R 3 f in i ϕ dpdx - R + η ′ R 3 R 3 f i ϕ dpdxdt - R + η R 3 R 3 f i A • ∇ (x,
R 3 R 3 f i (A•∇ (x,p) +λT )(ϕ-ϕ ) dpdx = R 3 R 3 f i (A • ∇ (x,p) + λT )(ϕ -ϕ ) dpdx = 0.
Finally the formulation (84) is equivalent to -η(0) Since div (x,p) A + T λ = 0 and T f i = 0 the previous equation reduces to [START_REF] Goudon | Hydrodynamic limit for the Vlasov-Navier-Stokes equations. I. Light particles regime[END_REF].

R 3 R 3 f in i ϕ dpdx - R + η ′ R 3 R 3 f i ϕ dpdxdt (85) - R + η R 3 R 3 f i (A • ∇ (x,

[ 6 ] 3 . 2

 632 Proposition Assume that inf x∈R 3 B(x) > 0. Then T restricted to ker • T : D(T ) ∩ ker • → ker • , is a one to one map onto ker • . Its inverse belongs to L(ker • , ker • ) and we have the following Poincaré inequality u ≤ 2π ω 0 T u , ∀u ∈ D(T ) ∩ ker •

  r, z) depending on time t and the invariants x, r = |p ∧ b(x)|, z = (p • b(x))

v 1 (

 1 cos(α + θ), ρ sin(α + θ), x 3 + θ -(α + θ), r, z)dα cos α, ρ sin α, x 3 + θα, r, z) dα and therefore v 1 ∈ ker T 1 . The next propostion is similar to Proposition 3.1.

  2π[ we can write by using polar coordinates (ρ, θ) with angles measured with respect to α cos(θα), ρ sin(θα), x 3 , r, z)χ(ρ, x 3 + θα, r, z)ρ dx 3 dθdρ = cos(θα), ρ sin(θα), x 3 + α, r, z)χ(ρ, x 3 + θ, r, z)ρ dx 3 dθdρ.

g 1 e 1 e

 11 Plugging the termp m • ∇ x f = ξ • ∇ x,r,z g e in (48) we obtain the equation

  and after performing the changeF e (t, x, v) = G e (t, x, r = |v ∧ b(x)|, z = v • b(x) ) the limit model (72) reduces to

)

  where µ(x, v) = m|v ∧ b(x)| 2 /2B(x), k(x) = |∂ x b b(x)| is the curvature of the magnetic lines and n(x) = ∂ x b b(x)/k(x) is the first normal to the magnetic lines.

  f εn i ω c (p ∧ b) • ∇ p ϕ dpdxdt = 0 implying that T f i (t) = 0, t ∈ R + .Therefore there is a functiong i = g i (t, x, r, z) such that f i (t, x, p) = g i (t, x, |p ∧ b(x)|, (p • b(x))).We use now the weak formulation of[START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] withC 1 test functions η(t)ϕ(x, p) with compact support in the open set {(t, x, p) ∈ R + × R 3 × R 3 : p ∧ b(x) = 0}and such that T ϕ = 0. Notice that these test functions can be written η(t)ψ(x, |p ∧ b(x)|, (p • b(x))), where ψ = ψ(x, r, z) are C 1 functions with compact support in R 3 × R ⋆ + × R. We obtain -η(0)

r 2 -∂ r 2πrg i zr 2m div x b +∂ z 2πrg i r 2

 22 2m div x b + e(b • E) ∂ z ψ dνdt = 0 saying that g i (0) = g in i and ∂ t (2πrg i )+div x 2πrg i zb m 2m div x b + e(b • E) = 0which is equivalent to (32). Since the magnetic momentum r 2 /(2mB(x)) is left invariant by (32), it is easily seen that (32) (supplemented by initial condition) is well posed in the phase-space (x, r, z) ∈ R 3 × R ⋆ + × R without any boundary condition at r = 0.In order to write a Vlasov equation for the dominant ion distribution f i in the phase space (x, p), we express the derivatives of ψ with respect to the derivatives of ϕ in (83)∂ z ψ = b(x)•∇ p ϕ, ∂ r ψ = p -(p • b)b |p ∧ b| •∇ p ϕ, ∇ x ψ = ∇ x ϕ-( ⊥ p•∇ p ϕ) t ∂ x b p |p ∧ b| , p∧b = 0leading to the following weak formulation, valid only for smooth test functions η(t)ϕ(x, p) vanishing near p ∧ b(x) = 0 and satisfying the constraint T ϕ = 0

R 3 R 3 AR 3 R 3 aR 3 R 3 aR 3 R 3 divR 3 R 3 (R 3 R 3 (R 3 R 3 χR 3 R 3 f

 3333333333333333 p) ϕ dpdxdt = 0(84) where A(t, x, p) = b(x) ⊗ b(x) p m , eb(x) ⊗ b(x)E(t, x) + ω i (x, p) ⊥ p .We intend to introduce a Lagrange multiplier in order to eliminate the constraint T ϕ = 0. By (30), (31) we know thata • ∇ (x,p) ϕ = A • ∇ (x,p) ϕ, a(t, x, p) = p m , eE(t, x)for any smooth test function η(t)ϕ(x, p) vanishing near p ∧ b(x) = 0 and satisfying the constraint T ϕ = 0. We claim that div (x,p) A = 0. Indeed for any functionϕ ∈ C 1 c (R 3 × R 3) such that T ϕ = 0 we haveR 3 R 3 div (x,p) A ϕ dpdx = -• ∇ (x,p) ϕ dpdx = -• ∇ (x,p) ϕ dpdx = -• ∇ (x,p) ϕ dpdx = (x,p) a ϕ dpdx = 0saying that div (x,p) A = 0. By Proposition 3.2 we deduce that there is a unique λ = λ(x, p) with zero average such that T λ = -div (x,p) A. Observe that A • ∇ (x,p) + λT leaves invariant the zero average functions. Indeed, for any smooth functions χ, ϕvanishing near p ∧ b(x) = 0, satisfying χ = 0, T ϕ = 0 we have A • ∇ (x,p) χ + λT χ)ϕ dpdx = div (x,p) (χA) + T (λχ))ϕ dpdx = -R 3 R 3 χA • ∇ (x,p) ϕ dpdx = -R 3 R 3 χ a • ∇ (x,p) ϕ dpdx =a • ∇ (x,p) ϕ dpdx = 0 implying that A • ∇ (x,p) χ + λT χ = 0. We establish now a weak formulation for f i , valid for any smooth function, let say ϕ, vanishing near p ∧ b(x) = 0. Observe that ϕ also vanishes near p ∧ b(x) = 0 and therefore (84) holds true with the test function ϕ . It is easily seen that R 3 R 3 f in i (ϕϕ ) dpdx = 0, i (t, x, p)(ϕϕ ) dpdx = 0, t ∈ R + and since ϕϕ vanishes near p ∧ b(x) = 0 we have

  p) + λT )ϕ dpdxdt = 0 for any smooth function vanishing near p ∧ b(x) = 0, saying that f i (0) = f in i and∂ t f i + div (x,p) (f i A) + T (λf i ) = 0 in D ′ ({p ∧ b(x) = 0}).

  x, p)ϕ(x, p) dpdx = R 3 ×R + ×R g i (t, x, r, z)ψ(x, r, z) dν, t ∈ R + .

Acknowledgements. This work was supported under the ANR EGYPT team. The authors address thanks for their warm hospitality and support.


Nikolaos Vlassopoulos
email: nikolaos.vlassopoulos@inria.fr

Nazim Fatès
email: nazim.fates@inria.fr

Clustering Behavior of a Bio-inspired Decentralized Aggregation Scheme

This note reviews a bio-inspired scheme for aggregating autonomous agents in the absence of global communication or coordination, a problem that is known as Decentralized Gathering. We present results on the clustering behavior of the agents, as we vary the main parameter that controls the agents' aggregation. Our observations show that there exist two phenomenologically different behaviors, characterized by two different evolutions of the number of clusters with time. We relate these different behaviors to the coupling of two factors: a change in the scale of the interaction range of the agents and a change in the significance of the local fluctuations in the model.

Introduction

Assume that a large number of autonomous and identical agents are scattered on a plane, and that there is no global authority to coordinate their actions nor any means of global communication. The problem of gathering these agents in a small area is known as the Decentralized Gathering. This problem is known to be difficult in general and is even be impossible to solve exactly in some continuous space frameworks [START_REF] Prencipe | Impossibility of gathering by a set of autonomous mobile robots[END_REF]).

One approach to solving the decentralized gathering problem consists of imitating the behavior of the amoebae species Dictyostelium discoideum [START_REF] Fatès | Solving the decentralised gathering problem with a reaction-diffusion-chemotaxis scheme -social amoebae as a source of inspiration[END_REF]). The main characteristic of this approach is the existence of an active environment that conveys simple messages among the agents, which are called virtual amoebae. The agents interact with the environment by either initiating the transmission of a message or by detecting the existence of messages in their local environment. These two types of interaction are the building components for a stigmergic behavior.

The virtual amoebae aggregation scheme has been shown to exhibit a rich dynamical behavior [START_REF] Fatès | Solving the decentralised gathering problem with a reaction-diffusion-chemotaxis scheme -social amoebae as a source of inspiration[END_REF]; Vlassopoulos and Fatès (2010)) and to be robust. In this note, we focus on a Cellular Automaton-based (CA) instance of the aggregation scheme, as it has been described in [START_REF] Fatès | Solving the decentralised gathering problem with a reaction-diffusion-chemotaxis scheme -social amoebae as a source of inspiration[END_REF], and present a qualitative description of the two contrasting, clustering behaviors that can be observed in the model. As we will describe, these behaviors result from a change of scale on the interactions among the agents, from short-ranged to long-ranged. Interestingly, the aggregation behavior persists despite this change of scale.

Virtual Amoebae Aggregation Scheme Active Environment

The existence of an active environment simplifies the agent behavior, by delegating parts of its complexity to the environment and allows for "self-sustained" messages, that can travel arbitrary large distances. The active environment is modeled with a two-dimensional Greenberg Hastings reaction-diffusion cellular automaton (GHCA, see [START_REF] Greenberg | Pattern formation and periodic structures in systems modeled by reaction-diffusion equations[END_REF]). The CA consists of an array of cells of dimensions L × L, a set of cell states, Σ, a set of transition rules for the states and, for each cell, a set of cells that constitute its neighborhood N c .

In the GHCA, Σ = {M, . . . , 0}, where M is called the excited state, M -1, . . . , 1 are called the refractory states and 0 is the neutral state. A cell becomes excited only if it is neutral and if at least one of its neighboring cells is excited. An excited cell will become refractory in the next time step and then decrease its state until it reaches the neutral state. The dynamics of the GHCA involve "waves" composed of wavefronts of excited cells followed by refractory cells that extend outwards from an excitation. Most importantly, when two reaction-diffusion wavefronts meet, they annihilate.

Agents

For simplicity, we consider agents as particles that can read the states of the cell on which they reside as well as the state of the neighboring cells. The virtual amoebae behavior is then summarized as follows: If the state of the cell where an agent resides is 0 (neutral), then, at each time step, the agent initiates ("fires") a reaction-diffusion process with probability λ, by setting the state of the cell to M . If the cell is neutral and an excited neighbor is detected, the agent moves towards the excited neighbor, choosing randomly if more than one are excited. Otherwise, if the cell is in a refractory state, do nothing. Here, λ, the firing rate, is the most important parameter of the aggregation model. In our study, each cell can hold at most two agents. Increasing the cell capacity affects mainly the spatial dimensions of the clusters and, for high values of λ, the aggregation time. Figure 1 shows the aggregation process for two different values of λ.

One may thus wonder what is causing the agents to aggregate in both cases, where we see a completely different quantitative behavior. A partial answer to this question is: the presence of fluctuations, both in terms of the density and in terms of emission of reaction diffusion waves.

Clustering Behavior

In a previous work [START_REF] Vlassopoulos | How fast can Virtual Amoebae aggregate?[END_REF]) we have shown that there exists an optimal value of λ such that the aggregation time is minimized. The two different clustering behaviors became apparent while studying the aggregation properties of the model (Fig. 1). In both cases, the agents, given a sufficient amount of time, will aggregate to a single cluster, but as we can observe, this is accomplished by exhibiting two completely different sequences of intermediate macroscopic configurations. In the second (bottom) figure, where λ is large enough, we observe that the agents form small clusters that progressively merge into bigger ones. This process continues until there are a few large clusters that persist for a relatively large amount of time before merging into a single one. However, in the first (top) figure we can observe that the agents "collapse" into a single cluster, without going through the formation of intermediate, persisting, clusters.

From our experiments so far, we have observed that this transition, i.e. from (A) the formation of competing clusters and progressive merging into a single one to (B) the "collapse" of the agents into a single cluster, and inversely, appears to be continuous. One important remark is that high values of λ favor small-range interactions among the agents, in the sense that the distance a wave will manage to travel, and consequently, the number of agents it will manage to interact with, before it is annihilated by the presence of other waves in the environment decreases as λ increases. Accordingly, small values of λ will allow a wave to travel larger distances and interact with more agents before it is annihilated, and therefore can be considered as larger-range interactions. To sum up, the aggregation process persists, in spite of the scale changes and the different macroscopic behavior that results from these changes.

Fluctuations as a Source of Order?

What is the "driving force" of aggregation in the different scales we described? The common denominators that "destabilize" persisting clusters and cause agents to collapse in both behaviors are the fluctuations, in terms of density and in terms of emissions. The density fluctuations exist even for very small values of λ, but in this case, where the interaction wavelengths are greater than the grid size, they are not significant and the agents aggregate into one cluster. For high values of λ they become important and are the main reason for the generation of the initial small local clusters that will subsequently merge, until only one cluster remains, but also one of the reasons that cause cluster to merge, since that out of two clusters with (sufficiently) different number of amoebae, we expect that the larger one will emit more waves, in average. The fluctuations in the emission times are the driving force that causes the clusters to merge, for both small and large values of λ. However, it is interesting that we observe the same effective behavior of the system in different scales. More precisely, the same "forces" that cause two amoebae to merge into a cluster, will cause two clusters to merge into a larger one and so on, until only one cluster exists. The merging behavior seems to be similar at different scales, which leaves us with the question: are there quantities that are invariant with respect to rescaling?

To conclude, we described a bio-inspired model that shows how it is possible to exploit the presence of fluctuations in a constructive way, in order to drive a system to a desired final state. The existence of an active environment simplifies the overall model design, but also increases the significance of fluctuations, that constitute a major factor to the operation and robustness of the model.

Figure 1 :

 1 Figure 1: Aggregation instances for different values of λ. Top: λ = 1 • 10 -5 Bottom: λ = 8 • 10 -2 . Agents are shown with green and the reaction diffusion wavefronts with orange. The initial number of agents is 400.