P Degond

S A Hirstoaga
email: sever.hirstoaga@math.univ-toulouse.fr

M-H Vignal
email: mhvignal@math.univ-toulouse.fr

The Vlasov model under large magnetic fields in the low-Mach number regime

This article is concerned with the kinetic modeling, by means of the Vlasov equation, of charged particles under the influence of a strong external electromagnetic field, i.e. when ε 2 , the dimensionless cyclotron period, tends to zero. This leads us to split the velocity variable in the Vlasov equation into fluid and random components. The latter is supposed to have a large magnitude of order 1/ε (which corresponds to the low Mach number regime). In the limit ε → 0, the resulting model is a hybrid model which couples a kinetic description of the microscopic random motion of the particles to a fluid description of the macroscopic behavior of the plasma. The microscopic model is a first-order partial differential system for the distribution function, which is averaged over the ultra-fast Larmor gyration and the fast parallel motion along the magnetic field lines. The perpendicular component (with respect to the magnetic field lines) of the bulk velocity is governed by the classical relations describing the E × B and diamagnetic drifts, while its parallel component satisfies an elliptic equation along the magnetic field lines.

Introduction

In this article, we are interested in the description of the dynamics of charged particles submitted to a large non uniform magnetic field. This problem is of great potential interest for describing strongly magnetized plasmas such as those encountered in Tokamak devices like ITER. The study of plasma confinement due to a large magnetic field requires to solve the Maxwell equations coupled to the description of the plasma turbulent transport. This transport can be modelled by using either a fluid description [START_REF] Beer | Toroidal Gyrofluid Equations for Simulations of Tokamak Turbulence[END_REF][START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF][START_REF] Dorland | Gyrofluid Turbulence Models with Kinetic Effects[END_REF][START_REF] Falchetto | Effect of Collisional Zonal-Flow Damping on Flux-Driven Turbulent Transport[END_REF][START_REF] Garbet | Global simulations of ion turbulence with magnetic shear reversal[END_REF][START_REF] Hammett | Developments in the Gyrofluid Approach to Tokamak Turbulence Simulations[END_REF][START_REF] Naulin | Electromagnetic transport components and sheared flows in drift-Alfvèn turbulence[END_REF][START_REF] Ottaviani | The gyro-radius scaling of ion thermal transport from global numerical simulations of ion temperature gradient driven turbulence[END_REF][START_REF] Scott | Free-energy conservation in local gyrofluid models[END_REF][START_REF] Xu | Low-to-high confinement transition simulations in divertor geometry[END_REF] or a kinetic description [START_REF] Bottino | Linear gyrokinetic simulations using particles for small perpendicular wavelength perturbations[END_REF][START_REF] Idomura | Global gyrokinetic simulation of ion temperature gradient driven turbulence in plasmas using a canonical Maxwellian distribution[END_REF][START_REF] Kim | Massively Parallel three dimensional Toroidal gyrokinetic Flux-Tube turbulence Simulations[END_REF][START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF][START_REF] Lin | Gyrokinetic simulations in general geometry and applications to collisional damping of zonal flows[END_REF][START_REF] Sugama | Comparison between kinetic and fluid simulations of slab ion temperature gradient driven turbulence[END_REF][START_REF] Tran | Global gyrokinetic simulation of Ion-Temperature-Gradient driven instabilities[END_REF]. Solving three dimensional fluid equations is certainly the less expensive way to solve the problem. However, a fluid description usually overestimates turbulent fluxes especially in a weakly collisional regime encountered in Tokamaks [START_REF] Dimits | Comparisons and physics basis of tokamak transport models and turbulence simulations[END_REF][START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. Indeed in Tokamaks, the plasma is carried to high temperatures. Thus, since the collision frequency decreases with increasing temperature, the plasma enters in a nearly collisionless regime.

By contrast, the kinetic model provides an appropriate description of turbulent transport in a fairly general context, but it requires to solve a six dimensional problem (3D in space and 3D in velocity) which leads to a huge computational cost. To reduce the cost of numerical simulations, it is classical to derive asymptotic models with a smaller number of variables than the kinetic description (see [START_REF] Ben Abdallah | On a hierarchy of macroscopic models for semi-conductors[END_REF] and references therein).

Here, we formally derive a new asymptotic model under both assumptions of large magnetic fields and low-Mach numbers. Large magnetic fields usually lead to the so-called driftkinetic limit (see [START_REF] Antonsen | Kinetic equations for low frequency instabilities in inhomogeneous plasmas[END_REF][START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF][START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF][START_REF] Hazeltine | Plasma Confinement[END_REF] for physics references and [START_REF] Bostan | The Vlasov-Poisson system with strong external magnetic field. Finite Larmor radius regime[END_REF][START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF][START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with a strong external magnetic field[END_REF][START_REF] Frénod | Long time behavior of the two-dimensional Vlasov equation with a strong external magnetic field[END_REF][START_REF] Frénod | Two-scale expansion of a singularly perturbed convection equation[END_REF][START_REF] Golse | The Vlasov-Poisson System with Strong Magnetic Field[END_REF] for mathematical results). In this regime, due to the large applied magnetic field, particles are confined along the magnetic field lines and their period of rotation around these lines (called the cyclotron period) becomes small. However, to our knowledge, the consideration of both large magnetic fields and low-Mach numbers is new. As we will see, considering low-Mach numbers brings a lot of interesting additional features.

We consider a simplified plasma model in which we focus on the dynamics of the ions. The coupling with the electrons is ignored and the electromagnetic field is assumed to be given. In future works, coupling the ion dynamics with those of the electrons and with the electromagnetic field is planned. To describe a collisionless ion dynamics (collisions can be neglected in Tokamaks in a first instance ; of course, there are situations where collisions must be included but we shall discard them in the present work), we use the Vlasov equation. In the large magnetic field regime, the Lorentz force term in the Vlasov equation is scaled by a large parameter, 1/ε 2 , where ε 2 stands for the dimensionless ion cyclotron period, i.e. the rotation period of the ion about a magnetic field line (or Larmor rotation). The so called drift-kinetic or gyro-kinetic regimes are reached when ε tends to zero (see [START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF][START_REF] Lifshitz | Physical Kinetics[END_REF]). We shall not dwell on the distinction between the drift and gyro kinetic regimes, as we are aiming at a different situation.

Indeed, in addition to the magnetic field being large, we assume that the ion mean velocity in the plasma is much smaller than the sound speed. In the Vlasov equation, this assumption implies that the fluid bulk velocity is much smaller than the magnitude of the random motion caused by thermal fluctuations [START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF]. Thus, we separate the slow scale of the fluid velocity from the fast scale of the random motion and we express the distribution function as a function of the random component of the velocity (which is also the velocity in the rest frame of the fluid) instead of the velocity in the laboratory frame. In the present work, we focus on the case where the Mach number is of order ε, i.e. it scales like the square root of the dimensionless cyclotron period. This scaling is natural since both quantities scale like the square root of the particle mass, and we can view both the large cyclotron freqency and small Mach numbers as a consequence of the small particle inertia. It also turns out that this scaling hypothesis gives rise to a rich structure in the asymptotic regime. The scaled Vlasov equation for the distribution of random velocities must be coupled to the fluid momentum equation which provides an equation for the bulk fluid velocity. Therefore, the unknowns of the Vlasov model in these new variables are the distribution of random velocities and the bulk fluid velocity. The goal of this paper is to investigate the limit ε → 0 of the Vlasov model in this scaling.

When ε → 0, the limit model consists of two sets of equations, one for the distribution function of random velocities, and one for the bulk fluid velocity. The distribution function of random velocities only depends on space, time and two components of the velocity, corresponding to the parallel component along the magnetic field line and the magnitude of the perpendicular velocity. In other words, the distribution function is independent of the (gyro)-phase of the perpendicular velocity in the plane normal to the magnetic field line. This is a consequence of the ultra-fast cyclotron rotation about the magnetic field lines. It is convenient to express the distribution of random velocities in terms of the parallel velocity and the magnetic moment (or adiabatic invariant), which is proportional to the perpendicular energy divided by the magnitude of the magnetic field. Now, the distribution function in these new variables satisfies a first order differential system with a constraint. A Lagrange multiplier allows to express this constraint in the differential system. The constraint expresses that the distribution function is constant along the trajectories of the fast parallel motion along the magnetic field lines. This motion is characterized by the constancy of the magnetic moment and of some kind of pseudo-energy in the parallel direction, which are the two adiabatic invariants of this motion. If a global change of variables from the phase space variables to the adiabatic invariants can be found, it is possible to eliminate the Lagrange multiplier of the constraint and to express the model as a transport equation in the space spanned by these invariants. This transport model describes how the electric and pressure forces as well as spatio-temporal variations of the magnetic field induce a slow evolution of the distribution function function in the space of adiabatic invariants. However, it is not always possible to find such a global change of variables, and, in this situation, the formulation of the problem as a constrained transport equation is the only possible expression of the system. Additionally, in most instances, it will provide a more flexible formulation for numerical discretization.

The equation for the bulk fluid velocity is split in two equations, one for the perpendicular component to the magnetic field lines, one for the parallel component. The perpendicular component is given through an identity which simply relates this velocity to the E × Bdrift and diamagnetic drift velocities. The equation for the parallel component is more unusual. It is an ellptic equation which expresses how, in the zero-Mach number limit, the parallel fluid velocity must adjust in order to guarantee that, at any time, the parallel components of the pressure and electric forces balance along the magnetic field lines. This elliptic equation is highly anisotropic because posed on each magnetic field line. It is obtained by obtained through expressing the constraint of zero parallel force thanks to the moments of the distribution function.

The derivation of the model roughly follows the following steps: we first proceed with formal expansions of the two unknowns in powers of the parameter ε (Hillbert expansion) and we keep the first three orders. The expansion of the momentum conservation equation readily leads to an equilibrium constraint expressing that, in the zero-Mach number limit, the pressure force must balance the electromagnetic force. This constraint provides an explicit relation for the perpendicular component of the bulk fluid velocity. The parallel component is only given implicitely through this equilibrium constraint. Finding an explicit equation for it requires some moments of the distribution function. The actual computation of this equation is slightly involved. Now, carrying the Hilbert expansion procedure through for the distribution function equation is best done if we change the random velocity variable into a coordinate system consisting of the parallel velocity, the energy, and the angle of rotation (or gyrophase) around the magnetic field line. Thanks to this coordinate change, we show that the leading order term of the distribution function does not depend on the gyrophase.

Next, we realize that, at each level of the expansion, we are led to inverting the gyrophase averaging operator [START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF][START_REF] Hazeltine | Plasma Confinement[END_REF]. We show that the inverse operator can only act on functions satisfying a specific solvability condition, namely that their gyrophase average is zero. We find the asymptotic model in abstract form by imposing this solvability condition successively to the various terms of the expansion, following the classical Hilbert expansion procedure of kinetic theory. Providing explicit expressions of the abstract operators appearing as outcomes of the Hilbert expansion procedures requires somehow tedious computations, most of which will be skipped and given for the reader's convenience in an appendix.

The remainder of the paper is organized as follows. In section 2, we present the scaling which expresses the assumptions of strong magnetic field and low Mach number regime. Then, we present and comment the main result of this article, namely the asymptotic model. In Section 3, by using Hilbert expansions we derive the equilibrium constraint for the leading order fluid velocity and we study the equations concerning the leading order distribution function. In Section 4, we write the abstract asymptotic model and we provide the main computational steps which lead to the explicit partial differential system for the limit distribution function. The conservative form of the model (main result) in terms of the magnetic moment is obtained in Section 4.3. Finally, in section 4.4, we give the explicit formula for the perpendicular fluid velocity and we obtain the elliptic equation for the parallel part.

2 The model, the scaling and the main result

The Vlasov equation in a strong magnetic field

We are interested in the dynamics of a single species of positively charged ions in the plasma. At this stage of the study, the coupling with the electrons is discarded and the electromagnetic field is supposed given. In future work, the model will be expanded by taking into account the coupling with the electrons and with a self-consistant electromagnetic field.

We are interested in finding the asymptotic limit of the Vlasov equation describing the dynamics of the ions when they are submitted to a large external magnetic field and where additionally, the thermal fluctuations of the velocity are large compared with the bulk fluid velocity.

Denoting by m the ion mass and by q the positive charge of the ion, we start from the Vlasov equation

∂ t f + v • ∇ x f + q m (E + v × B) • ∇ v f = 0, (2.1)
where

f ≡ f (x, v, t) is the distribution function and x ∈ Ω ⊂ R 3 x , v ∈ R 3
v , and t ∈ R + are respectively the position, velocity, and time variables. The position x is supposed to belong to an open domain Ω ⊂ R 3 . In addition, the electric field E ≡ E(x, t) and the magnetic field B ≡ B(x, t) are assumed to be given.

We supplement this equation with incoming boundary conditions on the boundary ∂Ω of Ω:

f (x, v, t) = f B (x, v, t), x ∈ ∂Ω, v • ν(x) < 0, (2.2)
where ν(x) is the outward unit normal to ∂Ω at x and f B is supposed given. f B represents the distribution function of incoming particles in the domain Ω. We also prescribe an initial datum

f (x, v, 0) = f I (x, v), x ∈ Ω, (2.3)
where f I is the distribution function of particles initially present inside the domain Ω. Next, we introduce a set of characteristic scales from which an appropriate scaling of equation (2.1) will be derived. Let x be a typical length scale of the problem and let v be the ion velocity scale (typically v is the thermal velocity of the ions, (2k B T /m) 1/2 , where k B is the Boltzmann constant and T is the temperature scale). The time scale is therefore t = x/v. We denote by B ∈ R + the characteristic magnitude of the applied magnetic field and by E = vB that of the electric field. Thus, we define the new variables and given fields by

x ′ = x/x, v ′ = v/v, t ′ = t/t, E ′ (x ′ , t ′ ) = E(x, t)/E, B ′ (x ′ , t ′ ) = B(x, t)/B.
Subsequently, letting f the distribution function scale, we introduce the new unknown

f ′ (x ′ , v ′ , t ′ ) = f (x, v, t)/f.
Inserting all these changes into (2.1) and dropping the primes for clarity, we obtain the dimensionless equation

∂ t f + v • ∇ x f + qBt m (E + v × B) • ∇ v f = 0. (2.4)
When the external magnetic field is assumed to be large, the rotation period of the ions about the magnetic field lines becomes small. Denoting by ω c = qB m the characteristic ion cyclotron frequency, we introduce the dimensionless cyclotron period

ε 2 = 1 t ω c .
Then, under this scaling, the Vlasov equation (2.4) for f = f ε takes the form:

∂ t f ε + v • ∇ x f ε + 1 ε 2 (E + v × B) • ∇ v f ε = 0, (2.5)
with intial and boundaray conditions still given by (2.2) and (2.3).

Splitting the drift and fluid velocities

In many instances, the ion mean velocity in the bulk plasma is much smaller than the sound speed [START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF]. To take into account this observation, we assume that the fluid ensemble velocity is much smaller that the random component of the particle velocity. Therefore, besides the magnetic field being large, another key assumption of the present work is that of a low Mach number. Because of the large magnetic field, the random component of the velocity undergoes a fast motion around the magnetic field lines. We will see that this drives the distribution function towards a state which is (at least on the average) isotropic in the plane normal to the magnetic field line. This contributes to a reduction of the dimension of the problem.

To implement this idea, we need to decompose the velocity into the fluid ensemble velocity and its kinetic part. For this purpose, we define the local density n ε (x, t) and fluid velocity u ε (x, t) as follows:

n ε (x, t) = R 3 f ε (x, v, t) dv, n ε u ε (x, t) = R 3 f ε (x, v, t) v dv.
(2.6)

The decomposition of the particle velocity is performed through the change of variables

v = u ε (x, t) + c/ε, (2.7)
where c/ε is the random component of the particle velocity. Here the random velocity has a large magnitude 1/ε which is expressed in the choice of the scaling. The distribution function expressed in these new variables is denoted by

f ε (x, v, t) = ε 3 g ε (x, c, t). (2.8)
The scaling factor ε 3 in front of g ε preserves the density:

g ε (x, c, t) dc = f ε (x, v, t) dv.
Because of the definition of the velocity u ε according to (2.6), g ε satisfies the constraint:

g ε (x, c, t) c dc = 0. (2.9)
Therefore, g ε does not carry any information about the evolution of the mean velocity itself. So, the transformed equation from (2.5) through the change of variables (2.7) will consist of two parts, an equation for g ε on the one hand and the momentum conservation equation which allows to determine u ε on the other hand.

After some easy algebraic manipulations, we find that (2.5) is equivalent to the following system:

∂ t g ε + u ε • ∇ x g ε -c • (∇ x u ε )∇ c g ε + 1 ε c • ∇ x g ε + ∇ x • P ε n ε • ∇ c g ε + + 1 ε 2 (c × B) • ∇ c g ε = 0, (2.10)
∂ t (n ε u ε ) + ∇ x • (n ε u ε ⊗ u ε ) + 1 ε 2 ∇ x • P ε -n ε (E + u ε × B) = 0, (2.11)
where ∇ x u ε stands for the Jacobian matrix (∂ x i u ε j ) i,j and P ε for the pressure tensor,

P ε = g ε c⊗c dc.
The model is supplemented with incoming boundary conditions:

g ε (x, c, t) = g ε B (x, c, t), x ∈ ∂Ω, (u ε + 1 ε c) • ν(x) < 0, (2.12)
where g ε B (x, c, t) is related to f B through (2.7), and with initial conditions (2.5). The only thing which remains to be proved is that u ε is the average velocity according to definition (2.6). This is a consequence that the constraint (2.9) is satisfied by (2.10), (2.11). To prove this, let us denote by I ε (x, t) = g ε (x, c, t) c dc. Then, multiplying (2.10) by c and integrating over c, we deduce

g ε (x, c, 0) = g ε I (x, c), x ∈ Ω, ( 2
ε 2 ∂ t I ε + (u ε • ∇ x )I ε + (I ε • ∇ x )u ε + (∇ x • u ε )I ε = I ε × B.
With the initial condition I ε (x, 0) = 0, we obtain (2.9). This shows the claimed equivalence.

Our goal is to study the asymptotic model formally obtained by taking the limit ε → 0 in (2.10), (2.11). We notice that taking the first moment in (2.10) and closing with P ε = p(n ε )Id, where n ε is the density of g ε , Id is the identity matrix and p(n) is a suitable isentropic pressure law, we find the isentropic Euler system with Lorentz force, under the scaling used in [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF]. This remark guided our choice of the present scaling of the Vlasov system. Indeed, it was shown in [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF] that this scaling allows to derive a sound drift-fluid model for magnetic plasma confinement.

Main result: the asymptotic model

We first introduce some notations. We assume that B does not vanish and we define the director of the magnetic field by b = B/|B|. For a particle with random component of the velocity c ∈ R 3 , we introduce its parallel component c with respect to the magnetic field and its magnetic moment µ, given by

c = c • b, µ = 1 |B| |c| 2 -c 2 2 .
We define also define the parallel and perpendicular components of the bulk velocity u, by

u = u b + u ⊥ , u = u • b, u ⊥ = b × (u × b) .
The quantity µ is the magnetic moment of the particle in its rotation motion about the magnetic field, i.e. the magnetic flux through the disk enclosed by the circular trajectory of the particle. Indeed, (|c| 2c 2 )/2 is the kinetic energy of the transverse motion to B and is proportional to the square of the transverse velocity. But the transverse velocity is equal to the Larmor radius r L times the gyrofrequency ω c and r 2 L is proportional to the surface S enclosed by the particle motion while ω c is proportional to |B|. Thus

|c| 2 -c 2 2 ∼ |c ⊥ | 2 ∼ r 2 L ω 2 c ∼ S |B| 2 ,
where c ⊥ = cc b and ∼ means proportionality. Thus, µ ∼ S |B| which is the magnitude of the flux of B through S, as announced.

For simplicity, we assume that the boundary ∂Ω of Ω is a magnetic surface (i.e. at any point on the surface, the magnetic field is tangent to it). As a consequence, a magnetic field line starting inside the interior of Ω does not intersect the boundary ∂Ω. We also assume that the closure Ω of Ω is a compact set. These assumptions are true for instance in most parts of a Tokamak device, since the magnetic field lines are either closed or dense on magnetic surface, depending whether the safety factor is rational or not. We recall that a Tokamak geometry is that of a torus and each magnetic surface at equilibrium is also a torus. The safety factor of a field line is the number of turns around the small circle of the torus which are necessary for one turn around the large one, following the field line. If the safety factor is rational, the field line is closed otherwise the field line is dense on a magnetic surface. The case of magnetic field lines intersecting the boundary can easily be dealt with but will be discarded for simplicity.

The main result of this paper, which will be proven in section 4 is the following:

Theorem 2.1 In the formal limit ε → 0, any solution (g ε , u ε ) of (2.10), (2.11) converges to (g, u) given as follows: first, there exists a function

G = G(x, µ, c , t) such that g(x, c, t) = G(x, µ, c , t).
Defining G = 2 π |B| G, the functions G, u ⊥ and u are solutions to the following problem: there exists a function

K = K(x, µ, c , t) such that ∂G ∂t + S † G + C † K = 0, (2.14)
C † G := ∇ x • (c G b) + ∂ ∂c ((B • Φ) G) = 0, (2.15)
u ⊥ = E × b |B| + 1 n|B| (b × ∇ x p ⊥ ) + p -p ⊥ n|B| f , (2.16
)

-3(b • ∇ x ) ∇ x • (p u b) +2 ∇ x • u ∇ x • (p b) b +(E • b) ∇ x • (nu b) + +(∇ x • b)∇ x • (-3p + p ⊥ )u b +p ⊥ (∇ x • b) 2 u = R 3 , (2.17)
where K = 2 π |B| k. The right-hand side R 3 , given by (4.37), does not depend on u . Finally, the operator S † is defined by

S † G := ∇ x • u + µ ∇ x × b -b × Φ + ( c 2 |B| -µ) f G + ∂ ∂c -(∇ x u) : (b ⊗ b) + µ∇ x • f + Φ • f c G + ∂ ∂µ -( ∂ ∂t + u • ∇ x ) ln |B| -∇ x • u + ∇ x u : (b ⊗ b) + 1 |B| ∇ x • (B × Φ) -( c 2 |B| ∇ x • f + Φ • f) µ G , (2.18)
with Φ, f and F given by

Φ = 1 |B| (F -µ∇ x |B|), f = b × (b • ∇ x )b, F = ∇ x • P n , (2.19)
and with the density n and the pressure tensor P defined by

n = G(µ, c ) dµ dc , P = p ⊥ (Id -b ⊗ b) + p b ⊗ b, (2.20
) Eq. (2.17) is an elliptic equation for the parallel component u of u which is posed along the magnetic field line. Indeed, only operators like ϕ → b • ∇ x ϕ or ϕ → ∇ x • (bϕ) are applied to u in (2.17). In the case of closed magnetic field lines, this problem can be uniquely solved. In the cae of dense field lines on a magnetic surface, it is a conjecture that this problem can be uniquely solved in the space of almost periodic functions. The study of this problem will be developped in future work.

p ⊥ = G(µ, c ) µ|B| dµ dc , p = G(µ, c ) c 2 dµ dc , (2.21
Both equations follow from the relation

- ∇ x • P n + E + u × B = 0. (2.23)
which expresses that, in the zero Mach limit, the total force must vanish. When resolving this equation in terms of u ⊥ , we immediatly get (2.16). However, the projection of this equation onto the direction parallel to b gives rise to an implicit constraint on u :

(- ∇ x • P n + E) • b = 0. (2.24)
To resolve this constraint into an explicit equation for u , one must take derivatives of this equation with respect to t and eliminate the time derivatives of n, p and p ⊥ which appear from this operation by using velocity moments of the distribution function equation (2.14). This computation will be detailed at section (4.4). Note that a similar procedure was used at the level of the drift-fluid limit in [START_REF] Degond | An Asymptotic Preserving Scheme for the Euler equations in a strong magnetic field[END_REF].

Trajectories of the fast parallel motion

To understand the constrained transport model (2.14), (2.15), it is useful to introduce the negative adjoint operator C to C † , given, for a function F = F (x, µ, c ) by:

CF := c b • ∇ x F + (B • Φ) ∂F ∂c . (2.25) Since ∇ x • B = 0, a simple computation gives C † G = C † (2π|B|G) = 2π|B|CG = |B|C( G |B| ). (2.26)
Therefore, the constraint (2.15) can be equivalently written

CG = 0, (2.27)
and expresses that G = G/(2π|B|) is constant on the characteristics of (2.15) or (2.26). These characteristics are all the curves defined parametrically by (x = X(τ

), µ = M(τ ), c = C (τ ))
where (X(τ ), M(τ ), C (τ )) satisfies the following ODE:

dX dτ = C b, (2.28
)

dC dτ = b • F -M∇ x |B| , (2.29
)

dM dτ = 0. (2.30)
Here, t is frozen and τ is a fictitious time which parametrizes the trajectory. The right-hand sides of (2.28)-(2.30) are evaluated at (X(τ ), M(τ ), C (τ )) and t.

The integral curves of the ODE system (2.28)-(2.30) define the leading order trajectories of the particles in their fast parallel motion along the magnetic field lines. Indeed, their parallel velocity is C as seen on (2.28). Similarly, (2.30) states that their magnetic moment is invariant in this motion. Finally, (2.29) shows that the parallel force is a combination of the projections onto the magnetic field lines of the electric force E on the one hand, and of the mirror force -µ∇ x |B| on the other hand (see [START_REF] Hazeltine | Plasma Confinement[END_REF]Chapter 4]). This last term takes into account the fact that, if the norm of the magnetic field increases along the particle motion, more energy is converted into the rotation motion about B, which decreases the magnitude of the parallel velocity. Symmetrically, if this norm decreases along the particle path, energy is transferred from the rotation to the translation motion and the parallel velocity increases. Now, the constraint (2.27) reflects the fact that the fast parallel motion along the magnetic field line is so fast that it instantaneously relaxes G to a constant along these trajectories.

More precisely, there are three time scales for a particle moving in a large magnetic field: the fastest time scale (which we will refer to as the 'ultra-fast' time scale) corresponds to the cyclotron or Larmor rotation period about the magnetic field. Then, the second fastest scale (referred to as the fast time scale) is the scale of the parallel motion along the magnetic field line, and which is described by system (2.28)-(2.30). Then, there is a slow time scale, which is that of the various drifts across the magnetic field lines, due to spatio-temporal variations of the magnetic field, or to the electric or pressure forces, etc. In this model, we focus on the slow time scale and assume that the ultrafast and fast time scales result in the homogenization (through time averaging) of the distribution function G. Therefore, as a consequence of the avearging over the ultra-fast time scale, G is independent of the gyrophase α. And because of the fast parallel motion, G is constant along the trajectories (2.28)-(2.30).

While expressing this averaging effect directly on the particle trajectories is difficult, the use of the kinetic model directly provides a way to do it by imposing constraints on the distribution function. This easier derivation reflects the fact that, to some extent, the distribution function describes the particle dynamics in a weak (or equivalently statistical) sense. Averaging the trajectories over some fast component is best done by looking at the evolution of an observable of the system which is constant over this fast motion. Saying that G is independent of the gyrophase and is constant along the trajectories of the fast parallel motion is equivalent to saying that there is an equal probability fo find particles with different gyrophases, or at different locations along these trajectories.

Invariants of the fast parallel motion

The fast parallel motion can be characterized by some invariants. Considering a particular magnetic field line, we introduce its curvilinear abscissa s and we write the function B • Φ = b • (F -µ∇ x |B|) locally along this trajectory as the derivative with respect to s of an effective potental V(s, µ): We immediately deduce two integrals of this motion:

B • Φ = - ∂V ∂s (s, µ). ( 2
M = Constant, C 2 2 + V(s, M) = Constant. (2.35)
The second relation expresses the conservation of a kind of parallel mechanical energy in the course of the fast parallel motion. Therefore, the phase-space trajectories of the fast parallel motion are supported by magnetic field lines for their space dependence, have constant magnetic moment and are such that their parallel velocity is related to the curvilinear abscissa along the magnetic field line and to the magnetic moment according to the second relation of (2.35). The magnetic moment is often referred to as the first adiabatic invariant while the second invariant in (2.35) is related to the second adiabatic invariant.

Indeed, the second adiabatic invariant is introduced when the fast parallel motion along the magnetic field line has a bounded spatial range. Suppose that we use the second relation of (2.35) to express C as a function of the curvilinear abscissa s and that this relation implies that C 2 can only be non-negative on a bounded interval [s 1 , s 2 ]. Then, the particle is bound to oscillate between the two mirror points s 1 and s 2 . This is because, at s 1 and s 2 , the total energy of the particle, (which is constant along the trajectory) is entirely in the form of perpendicular energy. Then, no energy is available for the parallel motion, and the only possibility for the particle is to bounce back to the regions where the perpendicular energy is lower. This periodic motion is referred to in the physics literature as the bounce motion. By classical mechanical calculations, the period of this motion, or bounce period, which is also called the second adiabatic invariant, is related to the expression appearing in (2.35). However, this expression has a general meaning, while the second adiabatic invariant is only meaningful in the case where the parallel motion is spatially bounded.

The Lagrange multiplier K

Now, K appearing in (2.14) is the Lagrange multiplier of the constraint (2.15). Indeed, let us consider a test function F (x, µ, c ) such that

CF = 0, (2.36)
i.e. such that F is constant along the trajectories of the fast motion. Multiplying (2.14) by F and integrating with respect to (x, µ, c ) leads to

d dt (G, F ) + (S † G, F ) = 0, (2.37)
where we have used the notation

(G, F ) = G F dx dµ dc = G F 2π|B| dx dµ dc . (2.38)
Here, Both G and F are constant along the fast trajectories, and it is equivalent to know G punctually or to know the collections of all (G, F ) for all test functions F satisfying (2.36). Indeed, the duality formula (2.38) defines a function BG:

2π(BG, F ) = (G, F ), ∀F such that CF = 0. (2.39)
Like G, the function BG is constant on the fast trajectories (note that G = 2π|B|G is not !). Now, by the same considerations, (S † G, F ) determines a function S † G. This operator is unambiguously defined by the relation:

2π(S † G, F ) = (S † G, F ), ∀F such that CF = 0. (2.40)
Is it difficult to write an explicit expression of the operators BG and S † G without a parametrization of the trajectories of the fast motion. As described in section (2.4.3), this requires a parametrization of the magnetic field lines, and a global definition of the two invariants (2.35). Such a global parametrization is not available in general. Section 2.4.5 below examines the particular case where such a global parametrization exists.

However, definitions (2.39) and (2.40) allow to eliminate the Lagrange multiplier (K) and to write (2.14) according to

∂ ∂t (BG) + S † G = 0, (2.41)
The introduction of the Lagrange multiplier is just a convenient way to express the contraint that G is constant along the trajectories of the fast motion, or equivalently, to remove the constraint (2.36) on the test function F for the weak form (2.37) of (2.41). Indeed, in weak form, system (2.14), (2.15) takes the form of a mixed formulation

d dt (G, F ) + (S † G, F ) + (C † K, F ) = 0, ∀F (2.42) (CG, H) = 0, ∀H (2.43)
where now, the two test functions F and H are unconstrained. This is the mixed type formulation of the constrained problem (2.37), (2.36). We note that K is defined up to a solution of C † K = 0. It does not seem obvious to single out a more pertinent choice of K among these possible solutions. However, this multiplicity does not affect the problem since whatever choice will lead to the same solution G.

Expression in the global coordinate system spanned by the invariants (if it exists)

To provide more explicit expressions of the operators BG and S † G, we assume that there exists a global change of variables

x → (y(x), s(x)), (2.44) such that y ∈ R 2 . A magnetic field line is defined by y(x) = Constant, and s(x) is the curvilinear abscissa along the magnetic field line. We denote by (µ, W), the two invariants defined by (2.35) and we also assume that W can be defined globally as a function W(x, µ, c ). When the vector-valued function Ψ(x, µ, c ) defined by

Ψ(x, µ, c ) = (y(x), µ, W(x, µ, c )), (2.45)
takes a fixed value (y 0 , µ 0 , W 0 ), the point (x, µ, c ) in phase-space spans a given trajectory of the fast motion denoted by C (y 0 ,µ 0 ,W 0 ) . We denote by dl (y 0 ,µ 0 ,W 0 ) the length element on this curve. Now, we obviously can write G as a function of (y, µ, W):

G(x, µ, c , t) = Ḡ(y(x), µ, W(x, µ, c )). ( 2

.46)

As a particular test function, we choose

F (x, µ, c , t) = δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )). ( 2

.47)

By the coarea formula [START_REF] Federer | Geometric Measure Theory[END_REF], the expression involving the delta distribution means that, for any smooth (for instance, continuous) test function φ, we have

φ δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc = C (y 0 ,µ 0 ,W 0 ) φ dl (y 0 ,µ 0 ,W 0 ) ∆ , (2.48)
where ∆ = det(dΨ dΨ * ), dΨ is the derivative of the above defined function Ψ, dΨ * is the matrix transpose of Ψ and the expression dΨ dΨ * means a matrix product. We note that for any function G like (2.46), the following identity holds true:

G δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc = Λ(y 0 , µ 0 , W 0 ) Ḡ(y 0 , µ 0 , W 0 ) (2.49)
where

Λ(y 0 , µ 0 , W 0 ) = δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc . (2.50)
is the coarea (colength) of C (y 0 ,µ 0 ,W 0 ) . The function F given by (2.47) obviously satisfies the constraint (2.36), and therefore, (2.40) leads to the following expression of S † G:

(Λ S † G)(y 0 , µ 0 , W 0 ) = S † G δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc (2.51)
Now, let us define by S x ∈ R 3 , S µ ∈ R, S ∈ R, the fields such that

S † G = ∇ x • (S x G) + ∂ ∂µ (S µ G) + ∂ ∂c (S G).
(2.52)

We write

S † G = ∇ (x,µ,c ) • (SG),
with S = (S x , S µ , s ). Then, bu using the Green formula and the chain rule:

(Λ S † G)(y 0 , µ 0 , W 0 ) = ∇ (x,µ,c ) • (SG)δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc = -G S • ∇ (x,µ,c ) δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc = -G∇ (y,µ,W) [δ((y, µ, W) -(y 0 , µ 0 , W 0 ))] | (y,µ,W)=Ψ(x,µ,c ) (∇ (x,µ,c ) Ψ • S) dx dµ dc = ∇ (y 0 ,µ 0 ,W 0 ) • G(∇ (x,µ,c ) Ψ • S)δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc (2.53)
For the last identity, we have used that,

∇ (y,µ,W) [δ((y, µ, W) -(y 0 , µ 0 , W 0 ))] = -∇ (y 0 ,µ 0 ,W 0 ) [δ((y, µ, W) -(y 0 , µ 0 , W 0 ))] .
Finally, we can define the vector fields

Σ y ∈ R 2 , Σ µ ∈ R, Σ W ∈ R and Σ = (Σ y , Σ µ , Σ W ), such that Σ y (y 0 , µ 0 , W 0 ) = 2π|B|(∇ x y • S x )δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc , (2.54)
Σ µ (y 0 , µ 0 , W 0 ) = 2π|B|S µ δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc , (2.55)
Σ W (y 0 , µ 0 , W 0 ) = 2π|B|(∇ (x,µ,c ) W • S)δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc , (2.56) and write:

S † G = Λ -1 ∇ (y,µ,W) • (Σ Ḡ). ( 2

.57)

Similar considerations lead to the fact that

B(x, µ, c , t) = B(y(x), µ, W(x, µ, c )), (2.58)
with

(Λ B)(y 0 , µ 0 , W 0 ) = 2π|B| δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dx dµ dc .
(2.59)

Finally, system (2.14), (2.14) is equivalent to the following system for Ḡ(y, µ, W, t):

∂ ∂t ( B Ḡ) + Λ -1 ∇ y,µ,W (Σ Ḡ) = 0. (2.60)
This system is a classical first order hyperbolic equation for Ḡ. Note that the existence of a global change of variables like (2.45) is not always guaranteed. Of course, a local one always exist (using the rectification theorem for vector fields for instance), but it may not be globally defined. Even if it is the case, the practical determination of this change of variables may be a challenging problem, which makes the use of these coordinates very often impossible. If the magnetic field evolves in time, this change of coordinates must be updated at any new time step, whch increases the computational difficulty of the problem even further.

Another comment is about the finiteness of the integrals (2.50), (2.54) to (2.56) and (2.59). If the magnetic field lines are closed, these integrals are on a compact set, and they are finite. If the magnetic field lines are not closed, but are dense on a magnetic surface, then, by imposing smoothness conditions on the test functions F , we deduce that the integrals are not defined on a particular magnetic field line, but rather on the whole magnetic surface.Then, we must change the definitions of the coefficients of the system accordingly. But with this change, the integrals become surface integrals on a compact manifold again (since we assumed from the very beginning that the closure Ω of the domain Ω was compact), and are finite as well. Of course, this requires to shift from one definition of the coefficients to another one and is very difficult to realize in practice. The weak mixed formulation is more suitable because it allows to automatically shift from line integrals to surface integrals when passing from closed magnetic field lines to dens ones.

Initial and boundary conditions

Let us first make some comments on the geometric assumptions of the boundary ∂Ω, namely that ∂Ω is magnetic surface. It implies that no magnetic field line originating from the interior of Ω intersects ∂Ω. If this were the case, then constaint (2.27) would impose that the value of G inside the domain along this field line would be specified by the boundary condition at the point where this magnetic field line interesects the boundary. This would result in the fact that the limit model would actually be posed on a subdomain Ω ′ ⊂ Ω consisting of the union of all magnetic field lines which do not intersect the boundary. In order to avoid any further complicated geometrical discussion, we wilil discard this situation and assume that ∂Ω is a magnetic field surface.

The model (2.14), (2.15) has to be supplemented with suitable boundary conditions on ∂Ω. The case where there exist global invariants (y, µ, W) to the trajectories of the fast motion will guide us in determining what are the proper boundary conditions. Indeed, in this case, the model reduces to the first order differential system (2.60) in the invariant space (y, µ, W).

We first specify some notations. The domain Ω×R + ×R where the variable (x, µ, c ) takes its values transforms into a domain O through the transformation to the (y, µ, W) variables (2.45). Its boundary is denoted by ∂O. The boundary conditions must be prescribed on the part of ∂O where the vector field Σ is incoming, i.e. Σ • λ < 0 where λ denotes the outward unit normal to ∂O at the considered point on the boundary.

Let us first assume that g ε B given by (2.12) is such that

g ε B → G B (x, µ, c , t), x ∈ ∂Ω, c • ν(x) < 0, (2.61)
as ε → 0. We also assume that G B satisfies CG B = 0, i.e. is constant alont the phasespace trajectories of the fast motion. This assumption is consistant since, by the geometric assumption on ∂Ω, magnetic field lines starting from a point on ∂Ω will remain on ∂Ω.

Then, we note the following identity:

(∇ (y,µ,W) • (Σ Ḡ), F ) O + ( Ḡ, Σ • ∇ (y,µ,W) F ) O = ∂O Ḡ F (Σ • λ) dΣ(y, µ, W), (2.62)
where ( F , Ḡ) O stands for the integral of F Ḡ over O and dΣ(y, µ, W) is the superficial measure on ∂O. Similarly, we have

(∇ (x,µ,c ) • (SG), F ) + (G, S • ∇ (x,µ,c ) F ) = ∂Ω×R + ×R GF (S x • ν) dS(x) dµ dc , (2.63)
where dS(x) is the superficial measure on ∂Ω and we recall the definition (2.38). By the calculation developed in section 2.4.5, the first and second terms of the left-hand side of (2.62) are respectively equal to the first and second terms of the left-hand side of (2.63), we deduce that the right-hand sides of these two formulas are equal. Splitting the boundary Γ = ∂Ω × R + × R into Γ + , Γ -, and Γ 0 , according to whether (S x • ν) is > 0, < 0 or = 0, and splitting ∂O into ∂O + , ∂O -, ∂O 0 in a similar way according to the sign of (Σ • λ), we get

∂O + Ḡ F (Σ • λ) dΣ(y, µ, W) - ∂O - Ḡ F |(Σ • λ)| dΣ(y, µ, W) = = Γ + GF (S x • ν) dS(x) dµ dc - Γ - GF |(S x • ν)| dS(x) dµ dc . (2.64)
We now make a very strong assumption, without which it is difficult to prescribe sound boundary conditions for this model. We assume that the sign of (S x • ν) is constant along the trajectories of the fast motion. Then, it is consistant to assume that the test function F satisfies the constraint (2.36), vanishes identically on Γ + but does not vanish and is actually arbitrary on Γ -. For such a test function, we deduce that

∂O - Ḡ F |(Σ • λ)| dΣ(y, µ, W) = Γ - GF |(S x • ν)| dS(x) dµ dc .
(2.65) Therefore, we prescribe the following boundary condition for the solution G of (2.14), (2.15):

Γ - GF |(S x • ν)| dS(x) dµ dc = Γ - 2π|B|G B F |(S x • ν)| dS(x) dµ dc .
(2.66) for any test function F which satisfies the constraint (2.36) and vanishes identically on Γ + , and where the datum G B is given by (2.61).

If it is possible to use global invariants as coordinates (y, µ, W), we can use (2.65) with a test function F of the form (2.47) with (y 0 , µ 0 , W 0 ) ∈ ∂O -. This leads to

( Ḡ|(Σ • λ)| dΣ)(y 0 , µ 0 , W 0 ) = = Γ - 2π|B|G B |(S x • ν)| δ(Ψ(x, µ, c ) -(y 0 , µ 0 , W 0 )) dS(x) dµ dc , , (2.67)
which gives an explicit prescription for Ḡ on the incoming boundary ∂O -.

The prescription (2.66), or in the cases where it is possible to use global invariants as coordinates, (2.67) are the prescribed boundary conditions for the model (2.14), (2.15).

The initial conditions are easier. We suppose that the initial conditions (2.13) are such that there exists G I (x, µ, c ), satisfying the constraint (2.27) such that

g ε I → G I , ∀(x, c) ∈ Ω × R 3 , (2.68)
as ε → 0 and we prescribe the boundary condition to be

G| t=0 = G I ., (2
.69)

Comparison with the literature

Let us introduce the following classical drift velocities due to the variations of the electromagnetic field [START_REF] Hazeltine | Plasma Confinement[END_REF]Chapter 4]:

V cd = c 2 || |B| f, the curvature drift due to the magnetic curvature (b • ∇ x )b, V gd = µ |B| (b × ∇ x |B|), the gradient drift, V ed = 1 |B| (E × b), the electric drift, V d|| = µ b • (∇ x × b) b, a drift parallel to the magnetic field.
With these definitions, the operator S † takes the following form

S † G := ∇ x • u || b + V d|| + V ed + V cd + V gd G + ∂ ∂c -b • ∇ x u -∇ x • (V gd + V d|| ) + (b • ∇ x )b • V ed + ∇ x |B| |B| • (V d|| -V gd ) c G + ∂ ∂µ - ∂ t |B| |B| -(b • ∇ x )b • V ed -∇ x • V ed + c 2 µ|B| ∇ x • V d|| - ∇ x |B| |B| • (V ed + V d|| ) µ G .
We recover the classical expressions of the spatial drift which can be found in the literature [START_REF] Hazeltine | Plasma Confinement[END_REF]Chapter 4]. By contrast, the drifts in c or µ which appear here do not appear in general in the literature. This work shows that, at least in the regime described by the proposed scaling, such drifts must be included otherwise, the limit model is not consistant with the Vlasov equation which we used as a starting point.

3 The asymptotic limit ε → 0: preliminaries

The Hilbert expansion

Our goal is to find the asymptotic limit ε → 0 of (2.10), (2.11). We start by assuming that g ε and u ε admit Hilbert expansions:

g ε = g 0 + εg 1 + ε 2 g 2 + . . . and u ε = u 0 + εu 1 + ε 2 u 2 + . . .
Inserting these expansions in equation (2.10), we find :

ε -2 term: -(c × B) • ∇ c g 0 = 0, (3.1) ε -1 term: -(c × B) • ∇ c g 1 = c • ∇ x g 0 + F 0 • ∇ c g 0 , (3.2)
ε 0 term: -(c × B) • ∇ c g 2 = ∂g 0 ∂t + u 0 • ∇ x g 0 -c • ((∇ x u 0 )∇ c g 0 ) + +c • ∇ x g 1 + F 0 • ∇ c g 1 + F 1 • ∇ c g 0 . (3.3) with F 0 = ∇ x • P 0 n 0 , F 1 = ∇ x • P n 1 .
Similarly, in (2.11), we obtain:

ε -2 term: F 0 = E + u 0 × B, ε -1 term: F 1 = u 1 × B,
We will see that the following orders in the expansion of eqs. (2.10) and (2.11) are not required.

Next, we introduce some formal notations. We denote by L, T , A, and D the following operators:

Lg = -(c × B) • ∇ c g, (3.4) T g = c • ∇ x g + F 0 • ∇ c g, (3.5
)

Ag = ∂g ∂t + u 0 • ∇ x g -c • ((∇ x u 0 )∇ c g), (3.6
)

Dg = F 1 • ∇ c g. (3.7)
Then, eqs. (3.1)-(3.3) take the form

Lg 0 = 0, (3.8) Lg 1 = T g 0 , (3.9) Lg 2 = Ag 0 + T g 1 + Dg 0 .
(3.10)

Since the leading order term in the expansion involves the operator L, which describes the effect of a circular motion around the magnetic field lines at infinite angular velocity, we now specifically examine the properties of this operator.

Study of L

We investigate the solutions of Lg = 0 and more generally, in a second step, of those of Lg = h for a given h. Since L operates on c only, the x and t variables are omitted. In what follows, we will not seek precise statements about the functional spaces. We assume that the functions are as regular as needed for the following statements to be correct. We denote by b = B/|B| the director of the magnetic field. We assume that it is always defined in the region of interest i.e. that B does not vanish, and that it is as smooth as needed.

We first introduce some notations. Let e 1 , e 2 be two vectors such that the set {e 1 , e 2 , b} forms a direct orthonormal basis of R 3 (i.e. b = e 1 × e 2 ). We denote by ( c

⊥ | = (c 2 1 + c 2 2 ) 1/2 = (2e -c 2 ) 1/2
where e = |c| 2 /2 is the energy. Now, we introduce a coordinate system for c which is derived from the cylindrical coordinates. Any c such that c ⊥ = 0 is uniquely defined by the triple (e, c , α) with (e, c ) ∈ D and α ∈ S 1 where

D = {(e, c ) | e ≥ 0, - √ 2e ≤ c ≤ √ 2e}, (3.11)
and S 1 is the one-dimensional torus R/(2πZ), such that

e = 1 2 (c 2 1 + c 2 2 + c 2 3 ), (3.12) c = c 3 , (3.13) cos α = c 1 /(c 2 1 + c 2 2 ) 1/2 , sin α = c 2 /(c 2 1 + c 2 2 ) 1/2 , (3.14)
or conversely

c 1 = (2e -c 2 ) 1/2 cos α, (3.15
)

c 2 = (2e -c 2 ) 1/2 sin α, (3.16) c 3 = c .
(3.17)

The angle α is the gyrophase while |c ⊥ | is proportional to the gyroradius of the particle in the magnetic field. Let g(c) = g(e, c , α) be the expression of g in this coordinate system where the function g(e, c , α) is defined on D × S 1 . Saying that α ∈ S 1 means that g is 2π-periodic with respect to α. Hence, Lg = 0 implies that ∂g/∂α = 0, i.e. g = g(e, c ), which proves (3.18).

In order to solve eq. Lg = h, we introduce the gyroaveraging operator Π defined for every function g(c) by

Πg (e, c ) = 1 2π R 3 g(c) δ |c| 2 2 -e δ(c • b -c ) dc = 1 2π S 1 g(e, c , α) dα, (3.20)
for all (e, c ) ∈ D. Πg is nothing but the mean value of g over the phase α.

Since our asymptotic model will be obtained by gyroaveraging the system (3.8)-(3.10), the following properties of Π will be useful. Their proof is easy and is left to the reader: Lemma 3.2 For any function g(x, c, t), we have:

Π ∂g ∂x i = ∂ ∂x i (Πg) + ∂b ∂x i • ∂ ∂c (Π(cg)) , (3.21)
Π ∂g ∂t = ∂ ∂t (Πg) + ∂b ∂t • ∂ ∂c (Π(cg)) , (3.22)
Π ∂g ∂c i = ∂ ∂e (Π(gc i )) + b i ∂ ∂c (Πg) , (3.23)
We deduce the following properties, which will be useful in the sequel:

Lemma 3.3
We have:

ΠL = 0, (3.24) f ∈ ker L ⇐⇒ f = Πf, (3.25)
and, for any function g,

ΠT g = ∇ x + F 0 ∂ ∂e • Π(cg) + (∇ x b) : ∂ ∂c Π(c ⊗ c g) + (F 0 • b) ∂ ∂c Πg, (3.26)
ΠAg = ∂ ∂t + u 0 • ∇ x Πg + ∂ ∂t + u 0 • ∇ x b • ∂ ∂c Π(c g) -∇ x u 0 : ∂ ∂e Π(c ⊗ c g) -((∇ x u 0 )b) • ∂ ∂c Π(c g) + (∇ x • u 0 )Πg, (3.27)
ΠDg = F 1 • ∂ ∂e Π(c g) + (F 1 • b) ∂ ∂c Πg, (3.28)
where : denotes the contracted product of two tensors. Concerning the moments of g, we have the Lemma 3.4 For any function g(x, c, t), we have:

n = R 3 g dc = D
Πg(e, c ) 2π de dc , Let g be a function lying in ker L for all (x, t), i.e. such that g = G(x, |c| 2 /2, c • b, t). We have:

n = D G(e, c ) 2π de dc , (3.29)
P = R 3 g c ⊗ c dc = p ⊥ (Id -b ⊗ b) + p b ⊗ b, (3.30)
p ⊥ = D G(e, c ) e - 1 2 c 2 2π de dc , p = D G(e, c ) c 2 2π de dc . (3.31)
In particular, we have

∇ x • P = ∇ x p ⊥ + b • ∇ x (p -p ⊥ ) + (p -p ⊥ )(∇ x • b) b + (p -p ⊥ )(b • ∇ x )b. (3.32)
Lemma 3.5 Let h be given. Then, equation

Lg = h, (3.33
where K is arbitrary. K can be uniquely determined if we impose to g to satisfy the cancellation condition Πg = 0.

(3.36)

The so-defined unique solution of (3.33) is denoted by g = L -1 h and the operator L -1 is called the pseudo-inverse of L. We have

L -1 h(e, c , α) = 1 |B| 2π 0 Γ(α, ϕ) h(e, c , ϕ) dϕ, (3.37)
where the Green kernel Γ(α, ϕ) is given by

Γ(α, ϕ) =        ϕ 2π , if 0 ≤ ϕ < α, ϕ 2π -1, if α ≤ ϕ ≤ 2π. (3.38)
L -1 h does not depend on the particular choice of the basis vectors (e 1 , e 2 ). If B is a smooth function of (x, t), L -1 g is a smooth function of (x, c, t) in the open set where B does not vanish.

Proof: Acting Π onto eq. (3.33) and using (3.24) shows that (3.34) is a necessary condition. Conversely, using (3. [START_REF] Hammett | Developments in the Gyrofluid Approach to Tokamak Turbulence Simulations[END_REF], we see that any solution of (3.33) is of the form (3.35). The only thing to prove is that this formula provides a periodic function of α (otherwise it is not possible to invert the change to cylindrical coordinates). This is true precisely if and only if h satisfies the solvability condition (3.34). Applying the cancellation condition (3.36) allows to uniquely determine K:

K(e, c ) = - 1 2π|B| 2π 0 h(e, c , ϕ) (2π -ϕ) dϕ .
Inserting this expression into (3.35), we find the expression of L -1 h given by (3.37), (3.38). We note that Γ can be extended by periodicity into a function of ϕ defined on S 1 . Now, if we add any constant to Γ, we do not change the result of (3.37) due to the solvability constraint (3.34). Precisely, it is an easy matter to see that, changing the origin e 1 from which we measure the phase α amounts to add a fixed constant to α, or to add a fixed constant to Γ. This shows that, despite the use of an origin for the phase α in (3.37), L -1 actually does not depend on this choice. As a result, L -1 h only depends on b and has the same regularity as b, i.e. it is smooth in the domain where B is smooth and non-zero.

Remark 3.6

The choice of the vectors e 1 , e 2 in the proofs above can be arbitrary. In particular, they are not required to form a smooth vector field. We note that this choice is equivalent to choose an origin for the phase α. By contrast, c and e do not depend on e 1 , e 2 . We have shown that the pseudo-inverse L -1 , despite the fact that its analytical expression (3.37) looks dependent on the choice of a particular origin for the phase, actually does not depend on it. This remark is important because generating a smooth orthogonal basis e 1 , e 2 in the plane normal to b is always locally possible but may be globally difficult. Additionally, such a basis is not always associated to an admissible coordinate system. Indeed, for this to be possible, commutation relations are required. This obviously is not the case in general. [START_REF] Bostan | The Vlasov-Poisson system with strong external magnetic field. Finite Larmor radius regime[END_REF] The asymptotic model: derivation

The asymptotic model in abstract form

In this section we derive an asymptotic model for the limit g 0 of g ε by formally passing to the limit ε → 0 in (2.10)- (2.11). This model will be deduced by solving the sequence of equations appearing in the Hilbert expansion (3.8)-(3.10). First, by a simple application of Lemma 3.1, eq. (3.8) is easily solved by: Proposition 4.1 There exists a function G(x, e, c , t) such that

g 0 (x, c, t) = G(x, |c| 2 /2 , c • b, t). (4.1)
Now, the goal is to find the equation satisfied by g 0 or G. For this purpose, we turn to (3.9). By a simple application of Lemma 3.5, we have: Proposition 4.2 Eq. (3.9) admits a solution g 1 if and only if g 0 satisfies the solvability condition ΠT g 0 = 0.

If this condition is satisfied, there exists a function k ∈ ker L (in other words, there exists a function K(x, e, c , t) and k = K(x, |c| 2 /2 , c • b, t)) such that

g 1 = L -1 T g 0 + k.
Note that, following (2.9) which has been proven equivalent to (2.11), we have:

g 1 c dc = 0.
Since, by construction, ΠL -1 = 0, and f dc = Πf dc for all functions f , we can write:

L -1 T g 0 c dc = ΠL -1 T g 0 c dc = 0.
Therefore, we deduce that

k c dc = 0.
Finally, again, by a simple application of Lemma 3.5, the equation satisfied by g 0 appears as the solvability condition of (3.10). For such a function, we compute Proposition 4.3 Eq. (3.10) admits a solution g 2 if and only if g 0 satisfies the equation

Π Ag 0 + T g 1 + Dg 0 = 0.
We now collect the model in the following theorem and discuss its properties. Theorem 4.4 The formal limit ε → 0 of problem (2.10), (2.11) leads to g ε → g 0 and u ε → u 0 , where g 0 and u 0 are solutions to the following abstract model: there exists a function k such that g0 = Πg 0 , k = Πk, (i.e. g0 and k are independent of α ),

(4.2) Π Ag 0 + T L -1 T g 0 + T k + Dg 0 = 0, (4.3) Π(T g 0 ) = 0, (4.4)
F 0 = ∇ x • P 0 n 0 = E + u 0 × B, F 1 = ∇ x • P n 1 = u 1 × B, (4.5)
k c dc = 0, (4.6)

where we recall that the tilde means that function is expressed in the coordinate system (3.12)- (3.14) and where the operators Π, A, T , L -1 and D are respectively defined by (3.20), (3.6), (3.5), (3.37) and (3.7).

Let us first note that eqs. (4.2) means that g0 and k do not depend of the gyrophase α. In this model, g 0 is determined by eq. ( 4.3) while the unknown function k plays the role of the Lagrange multiplier associated to the constraint (4.4). We will see in section 4.4 that u 0 is fully determined by the first equation of (4.5). Finally, it is not necessary to determine g 1 and u 1 . Indeed, g 1 appears in the definition of D but in the next section, using (4.5), we prove that ΠDg 0 = 0.

In the next sections, we make this model explicit in terms of a partial differential system.

∇ x + F 0 ∂ ∂e • (c Gb) + (∇ x • b) ∂ ∂c (e - 1 2 c 2 ) G + (F 0 • b) ∂G ∂c = 0.

The main equation (4.3)

We start with ΠAg 0 .

Lemma 4. [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] We have: We now denote γ 1 = L -1 T g 0 . With this definition, we have g 1 = γ 1 +k. We now compute ΠT γ 1 . By inspection of (3.26), we realize that we need to compute Π(c γ 1 ) and Π((c ⊗ c) γ 1 ) (we note that because of (3.36), Πγ 1 = 0). We first compute Π(c γ 1 ): Lemma 4.7 We have:

ΠAg 0 = ∂G ∂t + ∇ x • (u 0 G) -(∇ x • u 0 ) ∂ ∂e (e - 1 2 c 2 )G -b • ((∇ x u 0 )b) ∂ ∂e (-e + 3 2 c 2 )G + ∂ ∂c (c G) . ( 4
= ∂ ∂t + u 0 • ∇ x G + ∂ ∂t + u 0 • ∇ x b • ∂ ∂c (c Gb) -∇ x u 0 : ∂ ∂e G (e - 1 2 c 2 )(Id -b ⊗ b) + c 2 b ⊗ b -((∇ x u 0 )b) • ∂ ∂c (c Gb) + (∇ x • u 0 )G. ( 4
Π(c γ 1 ) = 1 |B| (e - 1 2 c 2 ) b × ∇ x + F 0 ∂ ∂e G + c ∂G ∂c (b • ∇ x )b . (4.11)
Proof: postponed in the Appendix.

Lemma 4.8 We have:

(∇ x b) : ∂ ∂c Π((c ⊗ c) γ 1 ) = = - 1 |B| ∂ ∂c (e - 1 2 c 2 )c b × ((b • ∇ x )b) • ∇ x + F 0 ∂ ∂e G . (4.

12)

Proof: postponed in the Appendix.

By collecting formulas (4.11) and (4.12) and inserting them into (3.26), we obtain the Lemma 4.9 We have:

ΠT γ 1 = ∇ x + F 0 ∂ ∂e • 1 |B| (e - 1 2 c 2 ) b × ∇ x + F 0 ∂ ∂e G + c ∂G ∂c (b • ∇ x )b - 1 |B| ∂ ∂c (e - 1 2 c 2 )c b × ((b • ∇ x )b) • ∇ x + F 0 ∂ ∂e G .
Since ∂ k/∂α = 0, the computation of ΠT k is the same as for ΠT g 0 : Lemma 4.10 We have:

ΠT k = ∇ x + F 0 ∂ ∂e • (c Kb) + (∇ x • b) ∂ ∂c (e - 1 2 c 2 ) K + (F 0 • b) ∂K ∂c .
Finally we calculate ΠDg 0 .

Lemma 4.11 We have: ΠDg 0 = 0.

Proof: Thanks to (3.7), (3.23) and (4.7), we get:

ΠDg 0 = F 1 • ∂ ∂e (Π(cg 0 )) + (F 1 • b) ∂ ∂c (Πg 0 ), = (F 1 • b) ∂ ∂e (c G) + ∂ ∂c G .
But, with (4.5), F 1 • b = 0. This concludes the proof.

The explicit form of the asymptotic model

From now on, for simplicity, we write u and F instead of u 0 and F 0 , respectively. We collect all the previous lemmas in the following theorem, which gives the explicit form of the drift kinetic model: Theorem 4.12 The formal limit ε → 0 of problem (2.10)- (2.11) leads to g ε → g 0 and u ε → u. There exists G such that g 0 (x, c, t) = g0 (x, e, c , α, t) = G(x, e, c , t), and G and u are solutions to the following partial differential system:

∂G ∂t + ∇ x • (uG) -(∇ x • u) ∂ ∂e (e - 1 2 c 2 )G -b • ((∇ x u)b) ∂ ∂e (-e + 3 2 c 2 )G + ∂ ∂c (c G) + ∇ x + F ∂ ∂e • 1 |B| (e - 1 2 c 2 ) b × ∇ x + F ∂ ∂e G + c ∂G ∂c (b • ∇ x )b - 1 |B| ∂ ∂c (e - 1 2 c 2 )c b × ((b • ∇ x )b) • ∇ x + F ∂ ∂e G + ∇ x + F ∂ ∂e • (c Kb) + (∇ x • b) ∂ ∂c (e - 1 2 c 2 ) K + (F • b) ∂K ∂c = 0, (4.13) ∇ x + F ∂ ∂e • (c Gb) + (∇ x • b) ∂ ∂c (e - 1 2 c 2 ) G + (F • b) ∂G ∂c = 0, (4.14)
F = ∇ x • P n = E + u × B, (4.15)
K c dc = 0, (

where K = K(x, e, c , t) is completely determined by the constraint (4.14) under some suitable boundary conditions, and where we recall the expressions of n and ∇ x • P are given in (3.29) and (3.32) respectively.

In the following section, we write the model in terms of more relevant variables.

The asymptotic model in the magnetic moment variable

In order to highlight the physical relevance of the model, it is useful to introduce the new variable (see [START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF])

µ = 1 |B| e - c 2 2 .
The quantity µ is the magnetic moment of the particle in its rotation motion about the magnetic field, see section 2.3. We introduce the change of variables G(x, e, c , t) = G(x, µ, c , t), K(x, e, c , t) = K(x, µ, c , t) .

We note that

(e, c ) ∈ D ⇐⇒ (µ, c ) ∈ D = R + × R, (4.17)
where D is defined by (3.11). Since B is divergence free, we have:

∇ x • b = B • ∇ x 1 |B| = -b • ∇ x (ln |B|).
Furthermore, the following formulas hold true:

∇ x G = ∇ x G -µ ∂G ∂µ ∇ x (ln |B|), ∂G ∂e = 1 |B| ∂G ∂µ , ∂G ∂c = ∂G ∂c - c |B| ∂G ∂µ , ∂G ∂t = ∂G ∂t -µ ∂G ∂µ ∂ ∂t (ln |B|).
From now on, we work in the variables (µ, c ) and we drop the overbars for clarity. From these formulas, we find Lemma 4.13 In the new variables, the constraint (4.14) is written

CG := c b • ∇ x G + b • F -µ∇ x |B| ∂G ∂c = 0. ( 4

.18)

A function G(x, µ, c , t) which satisfies the constraint (4.18) is constant along the curves where C is given by (4.18) and S by

(x = X(τ ), µ = M(τ ), c = C (τ ))
SG := u + µ ∇ x × b -b × Φ + ( c 2 |B| -µ) f • ∇ x G -(∇ x u) : (b ⊗ b) + µ∇ x • f + Φ • f c ∂G ∂c + - ∂ ∂t + u • ∇ x ln |B| -∇ x • u + ∇ x u : (b ⊗ b) + 1 |B| ∇ x • (B × Φ) - c 2 |B| ∇ x • f + Φ • f µ ∂G ∂µ , (4.21)
with Additionally, the velocity u satisfies the constraint:

Φ = 1 |B| (F -µ∇ x |B|), f = b × (b • ∇ x )b, F = ∇ x • P n ,
F = E + u × B. (4.22)
and K satisfies the constraint

K c dµ dc = 0, (4.23)
Proof: By performing the change of variables (4.17) into (4.13), we readily find that

SG := u • ∇ x G -∇ x u : (b ⊗ b) c ∂G ∂c - ∂ ∂t + u • ∇ x ln |B| + ∇ x • u -∇ x u : (b ⊗ b) µ ∂G ∂µ + ∇ x + Φ ∂ ∂µ • µb × ∇ x + Φ ∂ ∂µ G + ∇ x + Φ ∂ ∂µ • µc f ∂ ∂c - c |B| ∂ ∂µ G - ∂ ∂c - c |B| ∂ ∂µ µc f • ∇ x + Φ ∂ ∂µ G .
Now, we expand the second order operators and find:

∇ x + Φ ∂ ∂µ • µb × ∇ x + Φ ∂ ∂µ G = µ ∇ x × b -b × Φ • ∇ x G + 1 |B| ∇ x • (B × Φ) µ ∂G ∂µ , and
∇ x + Φ ∂ ∂µ • µc f ∂ ∂c - c |B| ∂ ∂µ G - ∂ ∂c - c |B| ∂ ∂µ µc f • ∇ x + Φ ∂ ∂µ G = c 2 |B| -µ f • ∇ x G + µ∇ x • f + Φ • f c ∂G ∂c - c 2 |B| ∇ x • f + Φ • f µ ∂G ∂µ ,
Lemma 4.16 The moment system satisfied by M m,q and K m,q is as follows:

∂ ∂t M m,q + ∇ x • u - b × F |B| M m,q + |B| (∇ x × b |B| ) -f M m,q+1 + f |B| M m+2,q + (m -q) (∇ x u) : (b ⊗ b) -(m -q) F |B| • f + q ∂ ∂t + u • ∇ x ln |B| + + q ∇ x • u - q |B| ∇ x • (b × F) M m,q + -m∇ x • f + (m -q) ∇ x |B| |B| • f + q |B| ∇ x • (b × ∇ x |B|) M m,q+1 + q |B| ∇ x • f M m+2,q +∇ x • (K m+1,q b) -m(b • F)K m-1,q + m(b • ∇ x |B|)K m-1,q+1 = 0, (4.26)
with the constraint

∇ x • (M m+1,q b) -m(b • F)M m-1,q + m(b • ∇ x |B|)M m-1,q+1 = 0, (4.27)
and with the convention that any moment with negative indexes is identically zero.

The infinite set of equations (4.26)-(4.27) for the positive integers m and q, contains the same information as the original model (2.14)- (2.15).

In particular, the first moments lead to Lemma 4.17 (i) The mass conservation equation holds:

∂n ∂t + ∇ x • (nu) = 0 . ( 4

.28)

(ii) We suppose that M 1,0 = 0 at the boundary of the domain. Then,

M 1,0 = G c dµ dc = 0 , (4.29)
everywhere.

(iii) The constraint (4.27) for m = 1 and q = 0 carries no information: it is redundant with the third equation in (2.19).

(iv) The pressures satisfy the following equations:

∂p ∂t + ∇ x • u - b × F |B| p + |B| (∇ x × b |B| ) -f M 2,1 + f |B| M 4,0 +2 (∇ x u) : (b ⊗ b) - F |B| • f p -2|B| ∇ x • f |B| M 2,1 +∇ x • (K 3,0 b) -2(b • F)K 1,0 + 2(b • ∇ x |B|)K 1,1 = 0, (4.30)
and

∂p ⊥ ∂t + ∇ x • u - b × F |B| p ⊥ + |B| (∇ x × b |B| ) -f |B|M 0,2 + f M 2,1 + -(∇ x u) : (b ⊗ b) + F |B| • f + ∇ x • u -∇ x • b × F |B| p ⊥ + ∇ x • f |B| |B|M 2,1 + |B|∇ x • (K 1,1 b) = 0. ( 4

.31)

Proof: (i) The equation for the density is obtained by letting m = q = 0 in (4.26):

∂ ∂t n + ∇ x • u - b × F |B| n + |B|(∇ x × b |B| ) -f p ⊥ |B| + f |B| p + +∇ x • (K 1,0 b) = 0. (4.32)
But, using (3.32), we get:

b |B| × nF = b |B| × ∇ x p ⊥ + f |B| (p -p ⊥ ).
Therefore, we have:

- b × nF |B| + |B|(∇ x × b |B| ) -f p ⊥ |B| + f |B| p = ∇ x × p ⊥ b |B| ,
and this term is canceled by the divergence operator in (4.32). With (2.22), eq. (4.28) follows.

Point (ii) follows from the application of the constraint (4.27) for m = 0 and q = 0. Indeed, we find

∇ x • (M 1,0 b) = 0 ,
out of which (4.29) follows from the assumption on the boundary conditions.

Point (iii) also follows from the inspection of the constraint (4.27) but with m = 1 and q = 0. Indeed, the left-hand side of this equation is

l.h.s. = ∇ x • (p b) -b • nF + b • ∇ x |B| |B| p ⊥ . ( 4
• nF = b • (∇ x • P) = b • ∇ x p + (p -p ⊥ )(∇ x • b) . Therefore, (4.33) is equal to l.h.s. = p ⊥ (∇ x • b) + b • ∇ x |B| |B| p ⊥ = p ⊥ |B| ∇ x • B = 0 ,
by the divergence free constraint on B. Therefore, the constraint (4.27) for m = 1 and q = 0 is redundant with the definition of F.

(iv) The equations for the pressures follow from the general moment equation (4.26) and (4.25). Now, we can turn towards the main result of this section, namely the Proof: In (4.24), we insert (4.34), the decomposition of the velocity into its aligned and transverse parts. By taking the vector product of (4.24) with b and using (3.32), we find that the transverse part of the velocity satisfies (2.16). The first component is the classical E × B drift. The two other components is the expression of the diamagnetic drift when the parallel and transverse pressures are different. We now turn to the difficult part: the determination of u . For it, we only have an implicit constraint, given by the projection of (4.24) onto b, i.e. Using the same methodology with eqs. (4.30), (4.31), we find

n(E

• b) -b • ∇ x p + (p -p ⊥ )(∇ x • b) = 0 . ( 4
∂p ∂t = -∇ x • (p u b) -2 p b • ∇ x u + R 1 , = -3∇ x • (p u b) + 2u ∇ x • (p b) + R 1 , ∂p ⊥ ∂t = -∇ x • (p ⊥ u b) + p ⊥ b • ∇ x u -p ⊥ ∇ x • (u b) + R 2 , = -∇ x • (p ⊥ u b) -p ⊥ u (∇ x • b) + R 2 ,
where

R 1 = -∇ x • p u ⊥ - b × F |B| + |B| (∇ x × b |B| -f ) M 2,1 + f |B| M 4,0 +2 -(∇ x u ⊥ ) : (b ⊗ b) + F |B| • f p + 2|B| ∇ x • f |B| M 2,1 -∇ x • (K 3,0 b) + 2(b • F)K 1,0 -2(b • ∇ x |B|)K 1,1 ,
and

R 2 = -∇ x • u ⊥ - b × F |B| p ⊥ + |B| (∇ x × b |B| ) -f |B|M 0,2 + f M 2,1 + (∇ x u ⊥ ) : (b ⊗ b) - F |B| • f -∇ x • u ⊥ + ∇ x • b × F |B| p ⊥ -∇ x • f |B| |B|M 2,1 -|B|∇ x • (K 1,1 b).
Inserting these formulas into (4.36), we find (2.17) with R 3 given by

R 3 = -(b • ∇ x )R 1 -(E • b)∇ x • (n u ⊥ ) -(∇ x • b) (R 1 -R 2 ) +n ∂ ∂t (E • b) - ∂b ∂t • ∇ x p -(p -p ⊥ )(∇ x • ∂b ∂t ). (4.37)
This is an elliptic equation for u which is invertible provided boundary conditions for u are given at the ends of the magnetic field line. which amounts to computing the m-th Fourier coefficients of h with respect to α. Then, the previous computation shows that

Π(c γ 1 ) = 1 |B|   -(2e -c 2 ) 1/2 Π 1 S T g 0 (2e -c 2 ) 1/2 Π 1 C T g 0 0   .
The third line corresponds to the applications of the cancellation condition (3.36).

We now need to explicitly compute T g 0 (so far, only ΠT g 0 was computed). Using the specific form (4.1) of g 0 , we find: Collecting these data, we deduce that:

Π 1 S T g 0 = 1 2 (2e -c 2 ) 1/2 (∇ x + F 0 ∂ e )G

Proof of Lemma 4.8

Using the same method as in the previous lemma, we show that in the basis (e 1 , e 2 , b), the matrix Π((c ⊗ c) γ 1 ) has the expression:

Π((c ⊗ c) γ 1 ) = 1 |4B|   -(2e -c 2 )Π 2 S T g 0 × × (2e -c 2 )Π 2 C T g 0 (2e -c 2 )Π 2 S T g 0 × -4(2e -c 2 ) 1/2 c Π 1 S T g 0 4(2e -c 2 ) 1/2 c Π 1 C T g 0 0   ,
where the symbol × indicates that the matrix is symmetric. So, again, we are left with the computation of Π 2 C,S T g 0 . In view of (5.2) and the fact that, obviously, Π 2 C,S (1) = Π 2 C,S (c) = 0, we need to compute Π 2 C,S (c ⊗ c). The same method as previously applies and leads to (5.6)

Now, we need to evaluate (∇ x b) : Π((c⊗c) γ 1 ). It is an easy matter to see that the contracted product of (∇ x b) with the first matrix of (5.6) is identically zero, as well as with the third one. The contracted product of (∇ x b) and the second matrix of (5.6) involves the expression

-(∇ x b) 31 (∇ x + F 0 ∂ e )G 2 + (∇ x b) 32 (∇ x + F 0 ∂ e )G 1 = = b • (∇ x + F 0 ∂ e )G × (b • ∇ x )b , = -(∇ x + F 0 ∂ e )G • b × (b • ∇ x )b .
(5.7)

Collecting (5.6) and (5.7) leads to the result (4.12) and ends the proof of the lemma.

)

  and the function K satisfies K c dµ dc = 0. (2.22) 2.4 Comments on the asymptotic model 2.4.1 The velocity u Eq. (2.16) provides an explicit relation for the transverse part of the velocity. In this expression, we recognize the E × B drift in the first term, and the diamagnetic drift in the second term. The third term is a drift term relating pressure anisotropy and the curvature of magnetic field lines. Indeed, (b • ∇ x )b is proportional to the curvature of the magnetic field line times the first normal to the curve. The vector b × (b • ∇ x )b is proportional to the curvature times the binormal to the curve (in the Frenet frame).

  1 , c 2 , c 3 ) the coordinates of c in this basis, with c 3 = c = c • b. We denote by c ⊥ = cc b the normal component of c to b. It has coordinates (c 1 , c 2 , 0) in this basis and its norm is given by |c

Lemma 3 . 1 c 2e -c 2 ,

 312 The null space ker L of L consists of functions which only depend on the parallel component c = c • b and on the energy e = |c| 2 /2, i.e. Lg = 0 ⇐⇒ ∃G(e, c ) with (e, c ) ∈ D, such that g(c) = G(|c| 2 /2 , c • b). and therefore, taking into account that |b × c| = |c ⊥ |, we obtain Lg = |B| ∂g ∂α . (3.19)

Proof: ( 3 .

 3 24) is immediately deduced from (3.23) and the fact that (B×c)•c = (B×c)•b = 0. (3.25) is obvious from the definition (3.20) and Lemma 3.1. The other formulas are simple applications of Lemma 3.2.

D

  Πg(e, c ) c 2π de dc .

4. 2

 2 The asymptotic model: explicit form 4.2.1 The constraint ΠT g 0 = 0 In this section, we consider the constraint (4.4): Proposition 4.5 Condition (4.4) is equivalent to the following equation for G:

Proof: 2 c 2 )

 22 We apply (3.26) with g = g 0 of the form (4.1). For such a function, an easy computation shows that Π(c g) = c Gb, (4.7) Π(c ⊗ c g) = G (e -1 (Idb ⊗ b) + c 2 b ⊗ b . (4.8) Then, we insert (4.7), (4.8) into (3.26) by noticing that (∇ x b) : (b ⊗ b) = ((b • ∇ x )b) • b = 0, since |b| = 1, and that (∇ x b) : Id = ∇ x • b.

. 10 )

 10 We first note that∂ ∂t + u 0 • ∇ x b • b = 0, since |b| = 1, which implies that the second term of (4.10) vanishes. Then, noting that (∇ x u 0 ) :Id = ∇ x • u 0 and (∇ x u 0 ) : (b ⊗ b) = ((∇ x u 0 )b)• b, the other terms combine into (4.9).

  and n and P being given by n = D G(µ, c ) 2π|B| dµ dc , P = p ⊥ (Idb ⊗ b) + p b ⊗ b, p ⊥ = D G(µ, c ) µ|B| 2π|B| dµ dc , p = D G(µ, c ) c 2 2π|B| dµ dc .

  .33) But, thanks to (3.30), (3.31), we get b

Lemma 4 .

 4 [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF] Let G, K and u satisfy eqs.(2.14),(2.15),(2.22) and (4.24) , thenu = u b + u ⊥ ,(4.34)with u = u • b and u ⊥ = b × (u × b) and u and u ⊥ are solutions to (2.16), (2.17).

  .35)To show how this leads to a well-posed equation for u , we take the time-derivative of (4.35) and use the continuity and pressure equations (4.28), (4.30), (4.31) to eliminate the time derivatives of n, p and p ⊥ . We first get from (4.35):∂n ∂t (E • b) + n ∂ ∂t (E • b)b • ∇ x ∂p ∂t + ∂b ∂t • ∇ x p + + ∂ ∂t (pp ⊥ ) (∇ x • b) + (pp ⊥ )(∇ x • p ⊥ ) (∇ x • b) = n ∂ ∂t (E • b) -∂b ∂t • ∇ x p -(pp ⊥ )(∇ x • ∂b ∂t ) . (4.36) Now, from (4.28), we have ∂n ∂t = -∇ x • (nu b) -∇ x • (nu ⊥ ) .

7 By ( 3 .-c 2 ) 1 / 2 1

 732121 [START_REF] Hazeltine | The drift kinetic equation for toroidal plasmas with large mass velocities[END_REF],(3.15)-(3.17) we have in the frame (e 1 , e 2 , b):Π(c γ 1 ) = 1 2π S 1 (cγ 1 )(e, c , α) dα, cos α (2ec 2 ) 1/2 sin α c   γ1 (α) dα. (5.1)But, thanks to (3.19), we get:1 2π S 1 cos α γ1 (α) dα = 1 2π S 1 d dα (sin α) γ1 (α) dα, cos α T g 0 (α) dα.We introduce the new averaging operators, for an arbitrary function h(c) and an arbitrary positive integer m: Π m S h(e, c ) = 1 2π S 1 sin(mα) h(e, c , α) dα, Π m C h(e, c ) = 1 2π S 1 cos(mα) h(e, c , α) dα,

T g 0

 0 = (F 0 • b) ∂G ∂c + c • ∇ x + F 0 ∂ ∂e G + (c ⊗ c) : ∇ x b T g 0 = Π 1 C,S (c) • ∇ x + F 0 ∂ ∂e G + Π 1 C,S (c ⊗ c) : ∇ x b ∂G ∂c ,and we are left with the task of computing Π 1 C,S (c) and Π 1 C,S (c ⊗ c). Using the same decomposition as for (5.1), we easily find: c) 1/2 c (e 2 ⊗ b + b ⊗ e 2 ), c ) 1/2 c (e 1 ⊗ b + b ⊗ e 1 ).

5 )

 5 b) 23 + (∇ x b) 32 , (5.3) c 2 ) 1/2 (∇ x + F 0 ∂ e )G 1 + c ∂G ∂c (∇ x b) 13 + (∇ x b) 31 , x + F 0 ∂ e )G) 2 ((∇ x + F 0 ∂ e )G) x b) 23 -(∇ x b) 32 (∇ x b) 13 + (∇ x b)From the fact that |b| = 1, we have (∇ x b)b = 0, out of which we deduce that (∇ x b) i3 = 0, for i = 1, 2, 3. The first vector in (5.5) can be easily identified with b × (∇ x + F 0 ∂ e )G while the second one, which reduces to (-(∇ x b) 32 , (∇ x b) 31 , 0) T (the exponent T denotes the transpose), is equal to b × ((b • ∇ x )b). Inserting these last remarks into (5.5) leads to(4.11) and ends the proof.

Π 2 Sc 2 4 (e 1 ⊗ e 2 + e 2 ⊗ e 1 ), Π 2 Cc 2 4 (e 1 ⊗ e 1 -e 2 ⊗ e 2 -c 2 4

 221222211224 (c ⊗ c) = 2e -(c ⊗ c) = 2ex b) 12 + (∇ x b) 21 ), Π 2 C T g 0 = Π 2 C (c ⊗ c) : ∇ x b ∂G ∂c = 2e ∂G ∂c ((∇ x b) 11 -(∇ x b) 22 ).With (5.3) and (5.4), we deduce that Π((c⊗c) γ 1 ) has the expression in the basis (e 1 , e 2 , b):Π((c ⊗ c) γ 1 ) = (2ec 2 ) x b) 12 -(∇ x b) 21 × × (∇ x b) 11 -(∇ x b) 22 (∇ x b) 12 + (∇ x b) ∇ x + F 0 ∂ e )G) 2 ((∇ x + F 0 ∂ e )G) x b) 32 (∇ x b) 31 0   .

Then, prescribing Ḡ on ∂Ω -is equivalent to prescribing the value of the integral at the lefthand side of (2.65) and, thanks to this relation, to prescribing the integral at the right-hand side of (2.65).

Acknowledgments. The authors wish to express their gratitude to G. Falchetto, X. Garbet and M. Ottaviani from the CEA-Cadarache and F. Deluzet from the Institut de Mathématiques de Toulouse, for fruitful discussions and encouragements. This work has been partially supported by the Marie Curie Actions of the European Commission in the frame of the 40 DEASE project (MEST-CT-2005-021122), by the CNRS and the Association Euratom-CEA in the framework of the contract 'Gyrostab' and by the CEA-Saclay in the framework of the contract 'Astre' # SAV 34160. This work was performed while the second author was an INRIA-Post-Doc at the Institut de Mathématiques de Toulouse.

where we have used that

and that

Collecting these various formulas together leads to (4.21).

We obtain the conservative form of the model given in Theorem 2.1 defining

Then

Additionally, K satisfies (2.22).

The relations between C † and S † on the hand and C and S on the other hand, are clarified in section 2.4.2.

Explicit equations for the limit fluid velocity

To find explicit equations for the velocity, we will need the moments of G. We introduce the general moments M m,q and K m,q of G and K, which are defined by (see [START_REF] Hazeltine | Plasma Confinement[END_REF]Chapter 6]) M m,q = G c m µ q dµ dc , K m,q = K c m µ q dµ dc .

We first note that n = M 0,0 , p = M 2,0 , p ⊥ = |B|M 0,1 .

(4.25)

Then, we obtain


Yi Zeng

Liyun Cao

Sheng Wan

Tong Guo

Shuowei An

Yan-Feng Wang
email: yswang@tju.edu.cn

Qiu- Jiao Du

Brice Vincent

Yue-Sheng Wang

Badreddine Assouar
email: badreddine.assouar@univ-lorraine.fr

Inertially amplified seismic metamaterial with an ultra-low-frequency bandgap

In last two decades, it has been theoretically and experimentally demonstrated that the seismic metamaterials are capable of isolating seismic surface waves. The inertial amplification mechanisms with small mass have been proposed to design metamaterials to isolate elastic waves in rods, beams and plates at low frequencies. In this letter, we propose an alternative type of seismic metamaterial providing an ultra-low-frequency bandgap induced by inertial amplification. A unique kind of inertially amplified

metamaterial is first conceived and designed. Its bandgap characteristics for flexural waves are then numerically and experimentally demonstrated. Finally, the embedded inertial amplification mechanism is introduced on a soil substrate to design a seismic metamaterial capable to strongly attenuate seismic surface waves around a frequency of 4 Hz. This work provides a promising alternative way to conceive seismic metamaterials to steer and control surface waves.

Seismic metamaterials (SMs), used to control seismic surface waves at low frequencies, have been investigated and proposed to protect critical infrastructures and ancient buildings in the last two decades. [START_REF] Meseguer | Rayleigh-wave attenuation by a semi-infinite two-dimensional elastic-band-gap crystal[END_REF][START_REF] Kim | Artificial seismic shadow zone by acoustic metamaterials[END_REF][START_REF] Brûlé | Experiments on seismic metamaterials: Molding surface waves[END_REF][START_REF] Miniaci | Large scale mechanical metamaterials as seismic shields[END_REF][START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF][START_REF] Zeng | Low-frequency broadband seismic metamaterial using I-shaped pillars in a half-space[END_REF][START_REF] Wang | Periodic in-filled pipes embedded in semi-infinite space as seismic metamaterials for filtering ultra-low-frequency surface waves[END_REF] The SMs composed of an array of artificial structures on or in a soil substrate usually possess the bandgap characteristics, whose working wavelengths is much longer than the size of the artificial structure. [START_REF] Miniaci | Large scale mechanical metamaterials as seismic shields[END_REF][START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF][START_REF] Zeng | Low-frequency broadband seismic metamaterial using I-shaped pillars in a half-space[END_REF][START_REF] Wang | Periodic in-filled pipes embedded in semi-infinite space as seismic metamaterials for filtering ultra-low-frequency surface waves[END_REF][START_REF] Aravantinos-Zafiris | Large scale phononic metamaterials for seismic isolation[END_REF][START_REF] Colombi | A seismic metamaterial: The resonant metawedge[END_REF][START_REF] Achaoui | Clamped seismic metamaterials: ultra-low frequency stop bands[END_REF] A number of experiments [START_REF] Meseguer | Rayleigh-wave attenuation by a semi-infinite two-dimensional elastic-band-gap crystal[END_REF][START_REF] Brûlé | Experiments on seismic metamaterials: Molding surface waves[END_REF][START_REF] Zeng | Broadband inverted T-shaped seismic metamaterial[END_REF] and simulations [START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF][START_REF] Du | Elastic metamaterial-based seismic shield for both Lamb and surface waves[END_REF][START_REF] Chen | Harnessing multi-layered soil to design seismic metamaterials with ultralow frequency band gaps[END_REF][START_REF] Liu | Ultra-wide band gap metasurfaces for controlling seismic surface waves[END_REF] have demonstrated that the SMs can significantly attenuate seismic waves at the frequencies of the bandgaps. In 1999, a phononic crystal in a marble was designed to experimentally demonstrate the existence of elastic bandgaps and the possible applications of SMs. [START_REF] Meseguer | Rayleigh-wave attenuation by a semi-infinite two-dimensional elastic-band-gap crystal[END_REF] Farhat et al. [START_REF] Farhat | Ultrabroadband elastic cloaking in thin plates[END_REF] proposed a design for a radially symmetric multilayered cloak to control flexural waves propagating in isotropic heterogeneous thin plates. The large-scale experiments, proposed by Brûlé et al., [START_REF] Brûlé | Experiments on seismic metamaterials: Molding surface waves[END_REF] have demonstrated that the SMs constituted of cylindrical holes periodically arranged in a soil substrate can attenuate seismic surface waves around a frequency of 50 Hz. It is easy to find evidence that the working frequency is not low enough to protect and shield human buildings. As Kadic et al. [START_REF] Kadic | Metamaterials beyond electromagnetism[END_REF] commented, the large-scale experiments are still far from being a functional cloak. Therefore, the SMs based on local resonance [START_REF] Liu | Locally resonant sonic materials[END_REF][START_REF] Sheng | Locally resonant sonic materials[END_REF][START_REF] Zeng | A broadband seismic metamaterial plate with simple structure and easy realization[END_REF][START_REF] Zeng | Subwavelength seismic metamaterial with an ultra-low frequency bandgap[END_REF] were proposed to obtain low-frequencies bandgaps for seismic surface waves. Finocchio et al. [START_REF] Finocchio | Seismic metamaterials based on isochronous mechanical oscillators[END_REF] introduced the mass-in-mass system in the soil to filter the harmonics of seismic waves at the fundamental resonance frequency of the building. Also, researchers proposed SMs composed of an array of pillars on a soil substrate to achieve low-frequencies bandgap, when they got inspiration that the seismic waves can be attenuated by the forests. [START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF][START_REF] Colombi | A seismic metamaterial: The resonant metawedge[END_REF][START_REF] Liu | Trees as large-scale natural metamaterials for lowfrequency vibration reduction[END_REF] To obtain wide bandgaps, different kinds of shape of the pillars [START_REF] Zeng | Low-frequency broadband seismic metamaterial using I-shaped pillars in a half-space[END_REF][START_REF] Zeng | Broadband inverted T-shaped seismic metamaterial[END_REF][START_REF] Lim | Built-up structural steel sections as seismic metamaterials for surface wave attenuation with low frequency wide bandgap in layered soil medium[END_REF][START_REF] Zeng | A Matryoshka-like seismic metamaterial with wide band-gap characteristics[END_REF] were proposed and the "rainbow trapping effects" [START_REF] Colombi | A seismic metamaterial: The resonant metawedge[END_REF][START_REF] Liu | Ultra-wide band gap metasurfaces for controlling seismic surface waves[END_REF][START_REF] Pu | Seismic metasurfaces on porous layered media: Surface resonators and fluid-solid interaction effects on the propagation of Rayleigh waves[END_REF][START_REF] Krödel | Wide band-gap seismic metastructures[END_REF][START_REF] Zeng | Coupling the first and second attenuation zones in seismic metasurface[END_REF] were used. However, at present, the bandgap for seismic surface waves below 5 Hz is still challenging to obtain just by using smallmass structures.

Inertial amplification mechanism, [START_REF] Yilmaz | Phononic band gaps induced by inertial amplification in periodic media[END_REF][START_REF] Yilmaz | Theory of phononic gaps induced by inertial amplification in finite structures[END_REF] which can achieve low-frequencies resonance just by using a small mass, has been introduced to design elastic metamaterials to obtain low-frequencies bandgaps. In 2006, Yilmaz et al. [START_REF] Yilmaz | Analysis and design of passive low-pass filter-type vibration isolators considering stiffness and mass limitations[END_REF] designed a passive low-pass filter-type vibration with low stop-band frequencies by using the lever-type structures. Then, the hybrid lever-type isolation systems with X-shaped structures were proposed to improve the band-stop characteristics at low frequencies. [START_REF] Liu | Vibration isolation using a hybrid lever-type isolation system with an Xshape supporting structure[END_REF][START_REF] Liu | Band stop vibration suppression using a passive X-shape structured lever-type isolation system[END_REF] In addition, the embedded amplification mechanism was proposed by Yilmaz et al. [START_REF] Yilmaz | Phononic band gaps induced by inertial amplification in periodic media[END_REF][START_REF] Yilmaz | Theory of phononic gaps induced by inertial amplification in finite structures[END_REF] to amplify the effective inertia of the wave propagation medium.

Frandsen et al. [START_REF] Frandsen | Inertial amplification of continuous structures: Large band gaps from small masses[END_REF] introduced the similar structures into a continuous elastic rod to obtain wide bandgap induced by inertial amplification. The three-dimensional metamaterials based on embedded inertial amplification mechanism have been realized.

It has been experimentally demonstrated that they are capable of isolating elastic waves in a wide frequency range. [START_REF] Taniker | Design, analysis and experimental investigation of three-dimensional structures with inertial amplification induced vibration stop bands[END_REF] Then, Yuksel et al. [START_REF] Yuksel | Realization of an ultrawide stop band in a 2-d elastic metamaterial with topologically optimized inertial amplification mechanisms[END_REF] have designed a two-dimensional solid structure with embedded inertial amplification mechanisms by using topology optimization to achieve an ultrawide stop band at low frequency. However, the embedded inertial amplification mechanisms have not been introduced into a half space to isolate surface waves.

In this letter, we propose a unique kind of SM with an ultra-low-frequency bandgap for seismic surface waves by using embedded inertial amplification mechanisms. At first, we design a kind of one-dimensional inertially amplified metamaterial (IAM). The unit cell of the IAM is composed of one base beam, two arms and three hinge joints. The characteristics of the bandgap induced by inertial amplification in the IAM is investigated by using experiments and simulations. Finally, the embedded inertial amplification mechanism is introduced to design the SMs for surface waves at ultralow frequencies. Ⅰ. [START_REF] Cao | Elastic Bound State in the Continuum with Perfect Mode Conversion[END_REF][START_REF] Zeng | Seismic metamaterials: Generating low-frequency bandgaps induced by inertial amplification[END_REF] The thickness t of the beam is 1 mm. The angle θ, which determines the magnitude of the amplification generated by the inertially amplified mechanisms, in this work is set at 5°. It is worth noting that hinge joint is not the necessary part for inertial amplification mechanisms. As shown in the literature, Yilmaz et al. [START_REF] Yilmaz | Analysis and design of passive low-pass filter-type vibration isolators considering stiffness and mass limitations[END_REF] designed a passive low-pass filter-type vibration with low stop-band frequencies by using the lever-type structures without hinge joint. However, the hinge joint is a simple way to achieve the inertial amplification mechanism in the metamaterials. Such as in the other literature, Frandsen et al. [START_REF] Frandsen | Inertial amplification of continuous structures: Large band gaps from small masses[END_REF] introduced the inertial amplification mechanism into a continuous elastic rod to obtain wide bandgap by using the hinge joints.

The simplified model of the inertially amplified mechanism of the unit cell of the IAM is illustrated in Fig. 1(d). The ma takes the value of the mass of one arm in the middle of the arm. The k and ka mean the stiffness of the beam and arm, respectively. Here, the ka can be regard as rigid at ultra-low frequencies. The angle θ between the arm and the beam is 5°. All the geometric parameters are shown in Figs. 1(c) and (e). It is noted that the lattice constant a is 55 mm and the unit of all the geometric parameters is millimeter.

In Fig. 1(d), the x1 and x2 represent the displacement of the input and output, respectively. At ultra-low frequencies, the horizontal and vertical displacement of the two ma are

𝑥 ! = (#$ ! %$ " ) ' ( 1
)
𝑧 ! = cot (𝜃) ($ ! ($ " ) ' (2) 𝑥 ) = ($ ! %#$ " ) ' (3)
𝑧 ) = cot (𝜃) ($ ! ($ " ) ' . (4)
Therefore, the kinetic energy (T) and potential energy (V) of the system are

𝑇 = * # +$̇# " %-̇# " . / + * # +$̇$ " %-̇$ " . / (5)
𝑉 = 0($ " ($ ! ) " / . ( 6
)
From the Lagrange's equation, i.e.,

0 = 1 12 / 34 3$̇" 0 - 34 3$ " , (7)
where the L is the Lagrange function:

𝐿 = 𝑇 -𝑉. (8)
Substituting Eq. ( 8) into Eq. ( 7), we obtain

/ /* # 56%789 " (:); <= 0 𝑥̈/ + 𝑘𝑥 / = /* # 5789 " (:)(/; <= 𝑥̈< + 𝑘𝑥 < . ( 9
)
When the two sides of Eq. ( 9) are equal to zero, the first resonance (ωr) and first antiresonance (ωar) [START_REF] Frandsen | Inertial amplification of continuous structures: Large band gaps from small masses[END_REF][START_REF] Yuksel | Realization of an ultrawide stop band in a 2-d elastic metamaterial with topologically optimized inertial amplification mechanisms[END_REF][START_REF] Li | Generating ultra wide low-frequency gap for transverse wave isolation via inertial amplification effects[END_REF] frequencies are obtained:

𝜔 > = 6 0 /* # [6%789 " (:)] <= ⁄ ( 10
)
𝜔 !> = 6 0 /* # [789 " (:)(/] <= ⁄ . (11)
As can be seen from Eqs. (10 -11), the ωr is always lower than the ωar. When the angle θ is less than 16.7°, the function 𝜂 = [5 + cot / (𝜃)] 16 ⁄ is larger than 1, which means the amplified motion for the mass is caused by the inertially amplified mechanism and the inertial mass is amplified η times. Of course, according to Eq. ( 10), if the stiffness k is weaker or the mass of two arms is heavier, the lower ωr can be obtained. It is worth noting that, comparing with the metamaterials based on local resonance, [START_REF] Mei | Dark acoustic metamaterials as super absorbers for low-frequency sound[END_REF][START_REF] Mei | Ultrathin metasurface with high absorptance for waterborne sound[END_REF] the inertial amplification mechanism is also based on the local resonance, which can be found from Eq. ( 10), but it can easily obtain lower-frequency bandgap because of the function 𝜂 = [5 + cot / (𝜃)] 16 ⁄ appeared in Eq. [START_REF] Achaoui | Clamped seismic metamaterials: ultra-low frequency stop bands[END_REF]. For example, in this work, the angle θ is 5°, so the inertial mass is amplified about 8.5 times. Therefore, the small mass can achieve strong resonance at low frequencies.

In the simulations, all results are calculated by using finite element software COMSOL 5.4 (Multibody Dynamics Module). The Floquet-Bloch boundary conditions are achieved on the unit cell along x direction by using both boundary similarity and pointwise constraint. In the Fig. 1(e), the hinge-connection condition is set on the hinge joints between two bodies. The transmission spectrums of the flexural waves in experiments and simulations are defined as

𝑇𝑆 = 20 × log <B (𝑅 ) 𝑅 ! ⁄ ) (12)
where Ra and Rb are the root mean square (RMS) of the out-of-plane velocities at the points a and b shown in Fig. 1(b), respectively.

𝜉 = E∫ 𝑤 𝑑𝐶 C E ∫ |𝑑𝑖𝑠𝑝. | 𝑑𝐶 C O ( 13
)
where w is the integral of the z component of the displacement fields in the yz-cut plane (to avoid the effect of the torsional modes in the beam); and disp. is the total displacement; and C means the whole unit cell. Therefore, in Fig. 2 Indeed, the experimentally measured attenuation of the IAM is far inferior compared to the simulated one. This is due to the fact that the hinge joints are rigid and smooth in the simulations, but not in the experiments. However, comparing with other two experimental samples (bare beam and TM), the attenuation zone is not difficult to find.

It is noted the unit cell of the TM has more mass than the IAM, but the frequency of the first bandgap of the TM for flexural waves is much higher than the IAM. The band structure of the TM, which is easy to calculate just by using Solid Module of the COMSOL, [START_REF] Brûlé | Experiments on seismic metamaterials: Molding surface waves[END_REF][START_REF] Miniaci | Large scale mechanical metamaterials as seismic shields[END_REF][START_REF] Zeng | Broadband inverted T-shaped seismic metamaterial[END_REF][START_REF] Ding | Metamaterial with simultaneously negative bulk modulus and mass density[END_REF] is not illustrated in this letter. At points B and D, almost all the displacement appears around the hinge joint between two arms, which means that there is a strong resonance only on the two arms. Based on this information, it is not difficult to infer that the low-frequencies bandgap is opened by the local resonance of the inertial amplification mechanism.

It is noted that the simplified model of embedded inertial amplification mechanism shown in Fig. 1(d) can also theoretically work on a half-space when the substrate provides the stiffness k. Therefore, in theory, the embedded inertial amplification mechanism is very suitable to design surface wave metamaterials, which, just with small mass, can induce low-frequency bandgap for surface waves. As shown in Fig. The lattice constant a in Fig. 3(a) is 2 m and the depth of the soil substrate is H = 500a [START_REF] Zeng | Broadband inverted T-shaped seismic metamaterial[END_REF][START_REF] Zeng | Seismic metamaterials: Generating low-frequency bandgaps induced by inertial amplification[END_REF]. The angle θ is also 5° and other geometer parameters are shown in Fig. 3(a).

The material parameters of the steel and soil used in this letter are shown in TABLE Ⅰ.

Figure 1 (

 1 Figure 1 (a) The samples used in the experiments (TM: traditional metamaterial; IAM: inertially

  Figure 1(c) illustrates the unit cell of the IAM considered in the numerical simulations

Figure 2 (

 2 Figure 2(a) illustrates the band structure of the IAM along x direction and the color bar

  (a), it is obvious that a bandgap induced by inertial amplification is generated for flexural waves in the frequency range from 70 to 108 Hz. As shown in Fig. 2(b), in the frequency range of the bandgap, a significant attenuation zone can be found in the simulation results. Comparing with all the experimental results, the IAM has an attenuation zone for the flexural waves around from 70 to 108 Hz.

Figure 2 (

 2 Figure 2 (a) The band structure of the IAM. The color bar represents the degree of the out-of-plane

Figure 2 (

 2 Figure 2(c) illustrates the vibration modes at the marked points A -D in the band

3

  (a), the embedded inertial amplification mechanism is introduced on a half-space to attenuate seismic surface waves. The unit cell of the seismic metamaterial (SM) is constituted of three parts: two steel arms represented by gray area, three hinge joints represented by dotted circle and soil substrate represented by yellow area. The blue and red lines represent Floquet-Bloch boundary condition and fixed constraint, respectively.

Figure 3 (

 3 Figure 3 (a) The unit cell of the SM with embedded inertial amplification mechanism. The blue and red

Figure 3 (

 3 Figure 3(b) illustrates the band structure of the one-dimension periodic SM. The sound

Figure 3 (

 3 Figure 3(d) illustrates the vibration modes at the marked points E -G in the band

TABLE Ⅰ :

 Ⅰ The material parameters used in this letter.

	Material	Density (kg/m 3 )	Young's modulus (Pa)	Poisson's ratio
	PLA	1086.3	3.4398 × 10 9	0.35
	Steel	7784	2.07 × 10 11	0.3
	Soil	1800	2 × 10 7	0.3

11991031, 11991032 and 41974059). The first author is grateful to the support of China Scholarship Council (Grant No. 202006250084).

DATA AVAILABILITY

The data that support the finding of the study are available from the corresponding author upon reasonable request.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.