Mihai Bostan
email: mbostan@univ-fcomte.fr

Gyrokinetic Vlasov equation in three dimensional setting. Second order approximation

Keywords: Vlasov equation, Guiding-center approximation, Average operator AMS classification: 35Q75, 78A35, 82D10

Gyrokinetic Vlasov equation in three dimensional

setting. Second order approximation

Introduction

Motivated by the confinement fusion, many research programs in plasma physics focus on strongly magnetized plasmas. It concerns the evolution of a population of charged particles under the action of strong magnetic fields B ε depending on some parameter ε > 0. Using the kinetic description and neglecting the collisions we are led to the Vlasov equation

∂ t f ε + p m • ∇ x f ε + q E(t, x) + p m ∧ B ε (x) • ∇ p f ε = 0, (t, x, p) ∈ R + × R 3 × R 3 (1)
with the initial condition

f ε (0, x, p) = f in (x, p), (x, p) ∈ R 3 × R 3 (2)
where f ε = f ε (t, x, p) ≥ 0 is the distribution function of the particles in the phase space (x, p) ∈ R 3 × R 3 , m is the particle mass and q is the particle charge. Generally we close the Vlasov equation by adding equations for the electro-magnetic field (E, B ε ) (i.e., the Maxwell equations or the Poisson equation). Here we consider only the linear problem (1), [START_REF] Bostan | The Vlasov-Poisson system with strong external magnetic field. Finite Larmor radius regime[END_REF] assuming that the magnetic field is stationary, divergence free and that the electric field derives from a given electric potential E(t) = -∇ x φ(t). We investigate the asymptotic behaviour of (1) when the magnetic field becomes large

B ε (x) = B(x) ε , B(x) = B(x)b(x), div x (Bb) = 0, 0 < ε << 1
for some scalar positive function B(x) and some field of unitary vectors b(x). We assume that B, b are smooth. Clearly the dynamics of the particles is dominated by the transport operator parallel to qB(x) p m ∧ b(x) • ∇ p . Assuming that ε is small enough we may expect that expansion like f ε = f + εf 1 + ε 2 f 2 + ... holds true and letting ε ց 0 it is easily seen that the leading order term belongs to the kernel of

T = qB(x) p m ∧ b(x) • ∇ p .
Notice that a family of independent invariants for T is given by x, |p ∧ b(x)|, p • b(x) and therefore the constraint T f = 0 is equivalent to

f (t, x, p) = g(t, x, r = |p ∧ b(x)|, z = p • b(x)).
Actually, plugging the above ansatz in [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations[END_REF] gives at the lowest order the divergence constraint T f = 0 and to the next order the evolution equation

∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f + T f 1 = 0. (3)
The key point is how to close [START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field. Guiding-center approximation[END_REF] with respect to the first order fluctuation density f 1 .

The idea is to project on the kernel of T by observing that the range of T is orthogonal to its kernel. Indeed, this will give a well-posed mathematical model, since we already know that f belongs to the kernel of T . The computations considerably simplify if we observe that the orthogonal projection on ker T is equivalent to averaging along the characteristic flow associated to T . The rigorous construction of the average operator (sometimes called by physicists the gyro-average operator in the context of gyrokinetic models) essentially relies on ergodic theory i.e., von Neumann's ergodic theorem [START_REF] Reed | Methods of Modern Mathematical Physics[END_REF] pp. 57. Employing this method we derive rigorously the guiding-center approximation in the three dimensional setting and we obtain the following Vlasov equation for the leading order particle density

∂ t f + b(x) ⊗ b(x) p m • ∇ x f + qb(x) ⊗ b(x)E + ω(x, p) ⊥ p • ∇ p f = 0 ( 4
)
where b(x) is the unitary vector field parallel to the magnetic field, the frequency ω(x, p) At the lowest order the particles are advected along the magnetic lines and only the parallel (with respect to b) electric field accelerates the particles. The plasma is confined along the magnetic lines and the transport operator in (1) becomes, in the limit ε ց 0

A x • ∇ x + A p • ∇ p = b(x) ⊗ b(x) p m • ∇ x + (qb(x) ⊗ b(x) E + ω(x, p) ⊥ p) • ∇ p .
But orthogonal drifts are expected at the next order. More general we investigate higher order approximations for [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations[END_REF], leading to a transport operator which takes into account the first order corrections

(A x + εA 1 x ) • ∇ x + (A p + εA 1 p ) • ∇ p .
Among these corrections we recover the electric cross field drift, the magnetic gradient drift and the magnetic curvature drift (cf. Theorem 5.2)

A 1 x = v ∧ + v GD + v CD + ... with v ∧ = E ∧ b B , v GD = |p ∧ b| 2 2m 2 ω c b ∧ ∇ x B B , v CD = (p • b) 2 m 2 ω c b ∧ ∂ x b b, ω c = qB m .
The main point is that the particle dynamics evolves on two time scales t and s = t/ε, the fast motion being associated to the large cyclotronic frequency 1 ε qB m . Accordingly the motion equations of the particles in (1) can be written

dX ε dt = P ε (t) m , dP ε dt = qE(t, X ε (t)) + 1 ε ω c (X ε (t))P ε (t) ∧ b(X ε (t)) (5)
where X ε (t) = X(t, t/ε) + εX 1 (t, t/ε) + ..., P ε (t) = P (t, t/ε) + εP 1 (t, t/ε) + ... . [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF] Plugging the above ansatz in ( 5) one gets at the lowest order ε -1

∂ s X = 0, ∂ s P = ω c (X)P ∧ b(X) (7)
and at the next order ε 0

∂ t X + ∂ s X 1 = P m (8)
∂ t P +∂ s P 1 = qE(t, X)+(∇ x ω c (X) •X 1 )P ∧b(X)+ω c (X) P ∧ ∂ x b (X)X 1 + P 1 ∧ b(X) .

From equation [START_REF] Brenier | Incompressible Euler and e-MHD as scaling limits of the Vlasov-Maxwell system[END_REF] we deduce that X, |P ∧ b(X)|, P • b(X) depend only on the slow time scale X = X(t), |P (t, s) ∧ b(X(t))| = R(t), P (t, s) • b(X(t)) = Z(t).

Moreover for any fixed t we have, by the second equation in [START_REF] Brenier | Incompressible Euler and e-MHD as scaling limits of the Vlasov-Maxwell system[END_REF] P (t, s) = cos(ω c (X(t)) s) b(X(t)) ∧ ( P (t, 0) ∧ b(X(t)) )

+ sin(ω c (X(t)) s) P (t, 0) ∧ b(X(t))

+ Z(t) b(X(t)) [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF] and therefore, at any fixed time t the momentum P is T c (X(t)) = 2π/|ω c (X(t))| periodic with respect to the fast variable s. Averaging the equation ( 8) with respect to s over one cyclotronic period T c (X(t)) one gets

dX dt = Z(t) m b(X(t)). (11)
At the leading order the particles are advected along the magnetic lines. Notice that

∂ s (P ∧ b(X)) = ∂ s P ∧ b(X) = ω c (X)(P ∧ b(X)) ∧ b(X) = ω c (X)( (P • b)b -P ) implying that ∂ s X 1 = P m -∂ t X = - ∂ s (P ∧ b(X)) ω c (X)m .
Therefore X 1 + P ∧b(X) mωc(X) is another invariant with respect to the fast motion. We can write

X(t) + εX 1 (t, s) = X(t) + ε X 1 + P ∧ b(X) mω c (X) -ε P ∧ b(X) mω c (X)
saying that during a cyclotronic period X(t) + εX 1 (t, t/ε) ≈ X ε (t) describes a circle of radius ε R(t) m|ωc(X(t))| into the plan orthogonal to b(X(t)). We compute now the acceleration along the magnetic lines by multiplying the equation ( 9) by b(X(t)) and averaging over one cyclotronic period (here • stands for the average with respect to s over one period).

For doing this observe that ∂ s P 1 • b(X(t)) = 0 and

∂ t P • b(X) = ∂ t (P • b(X)) -P • ∂ t b(X) = dZ dt - Z(t) m P • ∂ x b b(X) = dZ dt - Z 2 (t) m (b(X) • ∂ x b b(X)) = dZ dt . (12)
The average contribution of the electric force during a cyclotronic period is clearly

qE(t, X(t)) • b(X(t))
. It remains to compute the average contribution of the Laplace force. Notice that only the variation of the magnetic field direction accelerates the particles along the magnetic lines. Since X 1 + (P ∧ b(X))/(mω c (X)) is invariant with respect to the fast motion and div

x (Bb) = b • ∇ x B + B div x b = 0 we can write ω c (X)(P ∧ ∂ x b(X)X 1 ) • b(X) = 1 m ∂ x b(X) : (P ∧ b(X)) ⊗ (P ∧ b(X)) = |P ∧ b(X)| 2 2m (∂ x b : (I -b(X) ⊗ b(X))) = R 2 (t) 2m div x b = - R 2 (t) 2mB(X) ∇ x B • b(X).
We have obtained the diamagnetic force -µ(∇ x B • b) where µ(x, p) = |p∧b(x)| 2 2mB(x) is the magnetic moment. Combining the above computations yields

dZ dt = qE(t, X(t)) • b(X(t)) + R 2 (t) 2m div x b. (13)
Multiplying now the equation ( 9) by P/m and averaging over one cyclotronic period we deduce d dt

R 2 + Z 2 2m + ∂ s P 1 • P m = q( E(t, X(t)) • b(X(t)) ) Z(t) m + ω c (X) (P 1 ∧ b(X)) • P m .
Integrating by parts with respect to s one gets

∂ s P 1 • P = ω c (X) (P 1 ∧ b(X)) • P m .
Therefore the time variation of the cyclotronic momentum R(t) is given by

R(t) m dR dt + Z(t) m dZ dt = q(E(t, X(t)) • b(X(t))) Z(t) m .
Combining with [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field in quasineutral regime[END_REF] we obtain

dR dt = - Z(t)R(t) 2m div x b. ( 14
)
The dynamics of the particles with respect to the slow time scale in the phase space (x, r, z) is given by ( 11), ( 14), (13) leading to the limit model

∂ t g + z m b(x) • ∇ x g - zr 2m div x b ∂ r g + qE(t, x) • b(x) + r 2 2m div x b ∂ z g = 0 which is equivalent to (4) through the density change f (t, x, p) = g(t, x, |p ∧ b(x)|, p • b(x)).
The derivation of the second order approximation for (1) follows by employing similar techniques. Nevertheless it is a much difficult task, which requires complex computations, eventually the choice of appropriate coordinate system.

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

The guiding-center approximation for the Vlasov-Maxwell system was studied in [START_REF] Bostan | The Vlasov-Maxwell system with strong initial magnetic field. Guiding-center approximation[END_REF] by the modulated energy method, see also [START_REF] Bostan | High-electric-field limit for the Vlasov-Maxwell-Fokker-Planck system[END_REF], [START_REF] Brenier | Incompressible Euler and e-MHD as scaling limits of the Vlasov-Maxwell system[END_REF] for other results obtained by this method. The analysis of the Vlasov or Vlasov-Poisson equations with large external magnetic field have been carried out in [START_REF] Frénod | Homogenization of the Vlasov equation and of the Vlasov-Poisson system with strong external magnetic field[END_REF], [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field[END_REF], [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF], [START_REF] Frénod | The finite Larmor radius approximation[END_REF], [START_REF] Golse | The Vlasov-Poisson system with strong magnetic field in quasineutral regime[END_REF]. The numerical approximation of the gyrokinetic models has been performed in [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF] using semi-Lagrangian schemes. Other methods are based on the water bag representation of the distribution function: the full kinetic Vlasov equation is reduced to a set of hydrodynamic equations. This technique has been successefully applied to gyrokinetic models [START_REF] Morel | The water bag model and gyrokinetic applications[END_REF].

Our paper is organized as follows. In Section 2 we introduce the average operator and list its mathematical properties : orthogonal decomposition of L 2 functions into zero average functions and invariant functions along the characteristic flow, Poincaré inequality, etc. Section 3 is devoted to the derivation of the guiding-center approximation. This model is still a Vlasov equation. We investigare its conservative form and the geometry of its trajectories. We clearly identify invariants (magnetic moment, total energy) which allow us to reduce the dimension of the phase space. The asymptotic behaviour is studied in Section 4. We obtain both weak and strong convergence results. Section 5 is devoted to the second order approximation. One of the key points is to analyze the commutation properties between the average operator and first order differential operators.

Average operator

The main tool of our study is the average operator, which corresponds to the advection field dominating the transport operator in [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations[END_REF]. For simplicity we work in the

L 2 (R 3 × R 3
) framework but similar analysis can be carried out in any Lebesgue space

T u = div p (ω c (x) u p ∧ b(x)) , ω c (x) = qB(x) m D(T ) = {u(x, p) ∈ L 2 (R 3 × R 3 ) : div p (ω c (x) u p ∧ b(x)) ∈ L 2 (R 3 × R 3 )}.
We denote by • the standard norm of L 2 (R 3 × R 3 ). Notice that the above operator is local in x i.e., if u ∈ D(T ) then for a.a. x ∈ R 3 we have

u(x, •) ∈ L 2 (R 3 ) : div p (ω c (x) u(x, •) p ∧ b(x)) ∈ L 2 (R 3 ).
We denote by (X, P )(s; x, p) the characteristics associated to

ω c (x)(p ∧ b(x)) • ∇ p dX ds = 0, dP ds = ω c (X(s)) P (s) ∧ b(X(s)), (X, P )(0) = (x, p). (15)
Obviously X(s) = x for any s and by taking the scalar product of the second equation in [START_REF] Hazeltine | Plasma confinement[END_REF]

B(x) = ẽ1 (x) = b(x) ∧ (p ∧ b(x)) |b(x) ∧ p| , ẽ2 (x) = b(x) ∧ p |b(x) ∧ p| , ẽ3 (x) = b(x) .
Denoting ( P1 , P2 , P3 )(s) the coordinates of P (s) in the basis B(x) we obtain the equa-

tions d P1 ds = ω c (x) P2 (s), d P2 ds = -ω c (x) P1 (s), d P3 ds = 0
and therefore P3 (s) = P3 (0) Notice that the above formula holds also true in the case p ∧ b(x) = 0. The motions (X, P )(s; x, p) are T c (x) = 2π |ωc(x)| periodic for any initial condition (x, p) ∈ R 3 × R 3 . We introduce the average operator cf. [START_REF] Bostan | Transport equations with disparate advection fields[END_REF] u (x, p) = 1 T c (x)

P1 (s) = cos(ω c (x)s) P1 (0)+sin(ω c (x)s) P2 (0), P2 (s) = -sin(ω c (x)s) P1 (0)+cos(ω c (x)s) P2 (0). Taking into account the formula p = b(x) ∧ (p ∧ b(x)) + (b(x) • p) b(x) we deduce that P1 (0) = |p ∧ b(x)|, P2 (0) = 0, P3 (0)
Tc(x) 0 u(X(s; x, p), P (s; x, p)) ds = 1 2π S(x) u(x, |p ∧ b(x)| ω + (p • b(x)) b(x)) dω for any function u ∈ L 2 (R 3 × R 3 ), where S(x) = {ω ∈ S 2 : b(x) • ω = 0}.
Proposition 2.1 The average operator is linear continuous. Moreover it coincides with the orthogonal projection on the kernel of T i.e., u ∈ ker T :

R 3 R 3 (u -u )ϕ dpdx = 0, ∀ ϕ ∈ ker T .
Proof. For any function u ∈ L 2 (R 3 × R 3 ) we have for a.a.

x ∈ R 3 | u | 2 (x, p) ≤ 1 T c (x) Tc(x) 0 u 2 (x, P (s; x, p)) ds.
Taking into account that for any x ∈ R 3 the map p → P (s; x, p) is measure preserving one gets

R 3 R 3 u 2 (x, p) dpdx ≤ R 3 R 3 u 2 (x, p) dpdx saying that • ∈ L(L 2 (R 3 × R 3 ), L 2 (R 3 × R 3 )) and • L(L 2 (R 3 ×R 3 ),L 2 (R 3 ×R 3 )) ≤ 1.
It is well known that the kernel of T is given by the functions in L 2 invariant along the characteristics [START_REF] Hazeltine | Plasma confinement[END_REF]. Therefore we have

ker T = {u ∈ L 2 (R 3 × R 3 ) : ∃ v such that u(x, p) = v(x, |p ∧ b(x)|, (p • b(x)))}
Notice that for any u ∈ L 2 (R 3 ×R 3 ) its average u depends only on x, |p∧b(x)|, (p•b(x)).

Therefore u ∈ ker T . Pick a function ϕ ∈ ker T i.e.,

∃ ψ : ϕ(x, p) = ψ(x, |p ∧ b(x)|, (p • b(x))) ∈ L 2 (R 3 × R 3 )
and let us compute

I = R 3 R 3 (u -u )ϕ dpdx. Using cylindrical coordinates along b(x) axis yields I = R 3 R R + ψ(x, r, z) S(x) u(x, r ω + z b(x)) dω -2π u rdrdzdx = 0
and therefore u = Proj ker T u for any u ∈ L 2 (R 3 × R 3 ). In particular u = u for any u ∈ ker T and

• L(L 2 (R 3 ×R 3 ),L 2 (R 3 ×R 3 )) = 1.
For further use we inquire now about the solvability of T u = v. It is easily seen that if T u = v is solvable (i.e., v ∈ Range T ) then v = 0. Indeed, using the variational characterization of the average operator, we have for any function ϕ ∈ ker T

R 3 R 3 (v -0)ϕ dpdx = R 3 R 3 T u ϕ dpdx = - R 3 R 3 u T ϕ dpdx = 0
saying that v = 0. Generally we can prove that ker • = Range T . Indeed, since

• = Proj ker T and T ⋆ = -T we have ker • = (ker T ) ⊥ = (ker T ⋆ ) ⊥ = Range T .
Moreover we have the orthogonal decomposition of L 2 (R 3 ×R 3 ) into invariant functions along the characteristics [START_REF] Hazeltine | Plasma confinement[END_REF] and zero average functions

L 2 (R 3 × R 3 ) = ker T ⋆ ⊕ (ker T ⋆ ) ⊥ = ker T ⊕ Range T = ker T ⊕ ker • . ( 16
)
It happens that under additional hypotheses the range of T is closed, leading to the equality Range T = ker

u ∈ ker • ∩ D(T ) such that T u = v. For any α > 0 there is a unique u α ∈ D(T ) such that α u α + T u α = v. ( 18
)
Indeed it is easily seen that the solutions (u α ) α>0 are given by

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

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

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

u ∈ ker • ∩ D(T ) such that T u = v has compact support. Indeed, assume that supp v ⊂ {(x, p) ∈ R 3 × R 3 : |x| ≤ L x , |p| ≤ L p }. Since |P (s; x, p)| = |p| it is easily seen that supp u α ⊂ {(x, p) ∈ R 3 × R 3 : |x| ≤ L x , |p| ≤ L p }, α > 0 and therefore the weak limit u = lim n→+∞ u αn in L 2 (R 3 ×R 3 ) satisfies u ∈ ker • ∩D(T ), T u = v and supp u ⊂ {(x, p) ∈ R 3 × R 3 : |x| ≤ L x , |p| ≤ L p }.
Corollary 2.1 Under the hypotheses of Proposition 2.2 assume that the function v

belongs to W 1,∞ ([0, T ]; L 2 (R 3 × R 3 )) such that v(t) = 0 for any t ∈ [0, T ]. Then T -1 v belongs to W 1,∞ ([0, T ]; L 2 (R 3 × R 3
)) and we have

T -1 v W 1,∞ ([0,T ];L 2 (R 3 ×R 3 )) ≤ T 0 v W 1,∞ ([0,T ];L 2 (R 3 ×R 3 )) .
Proof. For any t ∈ [0, T ] we denote by u(t) the unique function of ker • ∩ D(T ) such that T u(t) = v(t). By Proposition 2.2 we deduce that

u L ∞ ([0,T ];L 2 (R 3 ×R 3 )) ≤ T 0 v L ∞ ([0,T ];L 2 (R 3 ×R 3 )) .
For any t ∈]0, T [ and h > 0 small enough we have

u(t + h) -u(t) ≤ T 0 v(t + h) -v(t) ≤ T 0 h v ′ L ∞ ([0,T ];L 2 (R 3 ×R 3 )) implying that u ∈ W 1,∞ ([0, T ]; L 2 (R 3 × R 3 )) and u ′ L ∞ ([0,T ];L 2 (R 3 ×R 3 )) ≤ T 0 v ′ L ∞ ([0,T ];L 2 (R 3 ×R 3 )) .
One of the crucial points when studying the asymptotic behaviour of ( 1) is how to propagate the regularity through the map T -1 . At the first sight this seems easily achieved by taking space/momentum derivatives in the equality T u = v and combining with the Poincaré inequality [START_REF] Littlejohn | Hamiltonian formulation of guiding center motion[END_REF]. Actually this arguments do not really work because the space/momentum derivatives may not commute with T and the average operator.

Indeed, notice that the Poincaré inequality provides an estimate for some derivative of u, let say ∂u, only if ∂u = 0. Therefore, if ∂ and • are not commuting, we may expect that ∂u = ∂ u = ∂0 = 0 and thus the Poincaré inequality can not be used. The above considerations lead naturally to derivatives along fields in involution with (0, ω c (x) p ∧ b(x)) ∈ R 6 i.e., fields c = c(x, p) ∈ R 6 such that the first order operator c(x, p) • ∇ x,p commutes with T . It was shown in [START_REF] Bostan | Transport equations with disparate advection fields[END_REF] that the average operator is commuting with derivatives along any such field in involution and finally we show that u inherits the regularity of v. It is also possible to appeal to a slightly different approach based on invariants [START_REF] Bostan | Transport equations with disparate advection fields[END_REF]. We recall that a complete family of invariants for

T is given by {x, |p ∧ b(x)|, p • b(x)}.
We will come back with more details about the propagation of regularity under the action of T -1 on zero average smooth functions, see Proposition 5.10.

Limit model

Using the properties of the average operator we investigate now the limit model of (1) when ε ց 0 by appealing to the method introduced in [START_REF] Bostan | Transport equations with disparate advection fields[END_REF] for general transport problems. We perform our computations by assuming high enough smoothness. We will see that the limit model is still a Vlasov equation, whose well posedness follows by standard arguments. We emphasize that the method we employ here has been studied in detail in [START_REF] Bostan | Transport equations with disparate advection fields[END_REF] (see also [START_REF] Bostan | The Vlasov-Poisson system with strong external magnetic field. Finite Larmor radius regime[END_REF]) for linear transport problems with even more general dominant advection field, with characteristic flow not necessarily periodic. The starting point consists in using a Hilbert expansion

f ε = f + εf 1 + ε 2 f 2 + ... (19)
Plugging the above ansatz into (1) yields

T f (t) = 0 (20) ∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f + T f 1 (t) = 0 (21) ∂ t f 1 + p m • ∇ x f 1 + qE(t, x) • ∇ p f 1 + T f 2 (t) = 0 (22)
. . .

From the constraint T f (t) = 0 we deduce that there is a function g = g(t, x, r, z) such that

f (t, x, p) = g(t, x, |p ∧ b(x)|, (p • b(x))). (23)
The time evolution of the dominant term f is described by (21) but we need to close this equation with respect to the first order correction f 1 . The key point here is to eliminate T f 1 by using the equality ker • = Range T . In this manner one gets the model

∂ t f + p m • ∇ x f + qE(t) • ∇ p f = 0. ( 24
)
Certainly we need to transform (24) into a more readable form taking into account the symmetries in (23). Notice that, by construction, the time derivative and the average operator are commuting. Therefore, since f = f ∈ ker T we obtain

∂ t f = ∂ t f = ∂ t f.
We compute now the averages of the derivatives with respect to space and momentum.

These computations become more complex due to the general geometry of the magnetic field. For the sake of the presentation we split them into separate lemmas.

Lemma 3.1 For any a, b ∈ R 3 , ω ∈ S 2 we have (ω • a) (ω • b) + (ω ∧ a) • (ω ∧ b) = a • b. In particular if a • b = 0 then (ω ∧ a) • (ω ∧ b) = -(ω • a) (ω • b). Proof. For any c ∈ R 3 we have |c| 2 = (ω • c) 2 + |ω ∧ c| 2 .

The conclusion follows

immediately by applying the previous formula with c ∈ {a, b, a + b}.

Lemma 3.2 Assume that f (x, p) = g(x, |p ∧ b(x)|, (p • b(x))). Then we have p m • ∇ x f = b(x)⊗b(x) p m •∇ x g- (p • b(x)) |p ∧ b(x)| 2m div x b ∂ r g+ |p ∧ b(x)| 2 2m div x b ∂ z g.
Proof. For any i ∈ {1, 2, 3} we have

∂ x i f = ∂ x i g + ∂ r g p ∧ b(x) |p ∧ b(x)| • (p ∧ ∂ x i b) + ∂ z g (p • ∂ x i b). Since |b(x)| = 1 we have ∂ x i b • b(x) = 0 and therefore by Lemma 3.1 one gets (b(x) ∧ p) • (∂ x i b ∧ p) = -(b(x) • p) (∂ x i b • p)
implying that

∂ x i f = ∂ x i g -∂ r g (p • b(x)) |p ∧ b(x)| (p • ∂ x i b) + ∂ z g (p • ∂ x i b).
We obtain the formula

p • ∇ x f = p • ∇ x g -∂ r g (p • b(x)) |p ∧ b(x)| ∂b ∂x : p ⊗ p + ∂ z g ∂b ∂x : p ⊗ p . Notice that (∇ x g)(x, |p ∧ b(x)|, (p • b(x))) ∈ ker T and therefore p • ∇ x g = p • ∇ x g = (b(x) • p) (b(x) • ∇ x g).
By direct computation we obtain that

p ⊗ p (x, p) = 1 2 |p ∧ b(x)| 2 (I -b(x) ⊗ b(x)) + (p • b(x)) 2 b(x) ⊗ b(x). (25)
Taking into account that ∂b ∂x : b

(x) ⊗ b(x) = b(x) • ∇ x |b(x)| 2 2 = 0 we deduce ∂b ∂x : p ⊗ p = 1 2 |p ∧ b(x)| 2 div x b
and finally one gets

p m • ∇ x f = b(x)⊗b(x) p m •∇ x g- (p • b(x)) |p ∧ b(x)| 2m div x b ∂ r g+ |p ∧ b(x)| 2 2m div x b ∂ z g. Lemma 3.3 Assume that f = g(x, |p ∧ b(x)|, (p • b(x))
) and E = E(x). Then we have

qE • ∇ p f = q(b(x) • E(x)) ∂ z g.
Proof. For any i ∈ {1, 2, 3} we have

∂ p i f = ∂ r g p ∧ b(x) |p ∧ b(x)| • (e i ∧ b(x)) + ∂ z g b i (x)
where B = {e 1 , e 2 , e 3 } is the canonical basis of R 3 . By Lemma 3.1 we know that

(b(x) • p) (b(x) • e i ) + |b(x) ∧ p| • (b(x) ∧ e i ) = p i and therefore |b(x) ∧ p| • (b(x) ∧ e i ) = (p -(b(x) • p) b(x)) i .
We obtain the formula

∇ p f = ∂ r g |p ∧ b(x)| (I -b(x) ⊗ b(x))p + ∂ z g b(x)
and finally, since

(∂ r g, ∂ z g)(x, |p ∧ b(x)|, (p • b(x))), E(x) ∈ ker T , one gets qE • ∇ p f = q(b(x) • E(x)) ∂ z g.
Combining the computations in Lemmas 3.2, 3.3 yields the following Vlasov equation in the phase space (x, r, z

) ∈ R 3 × R + × R ∂ t g + z m b(x) • ∇ x g - zr 2m div x b ∂ r g + r 2 2m div x b + q(b(x) • E(t, x)) ∂ z g = 0 (26)
whose characteristics (X, R, Z)(s; t, x, r, z) are given by

dX ds = Z(s) m b(X(s)) (27) dR ds = - Z(s)R(s) 2m div x b (X(s)) (28) dZ ds = R(s) 2 2m div x b (X(s)) + q(b(X(s)) • E(s, X(s))) (29) (X, R, Z)(t) = (x, r, z). ( 30
)
Certainly it is possible to write a Vlasov equation for the dominant distribution f . For this it is sufficient to express the derivatives of g with respect to the derivatives of f .

Recall that we have already obtained the formula

∂ t f = ∂ t g ∇ x f = ∇ x g-∂ r g (p • b(x)) |p ∧ b(x)| t ∂ x b p+∂ z g t ∂ x b p, ∇ p f = ∂ r g b(x) ∧ (p ∧ b(x)) |p ∧ b(x)| +∂ z g b(x).
Accordingly the derivatives of g write

∂ t g = ∂ t f, ∂ z g = b(x) • ∇ p f, ∂ r g = b(x) ∧ (p ∧ b(x)) |p ∧ b(x)| • ∇ p f ∇ x g = ∇ x f + b(x) ∧ (p ∧ b(x)) |p ∧ b(x)| • ∇ p f (p • b(x)) |p ∧ b(x)| t ∂ x b p -(b(x) • ∇ p f ) t ∂ x b p
leading to the following Vlasov equation

∂ t f + b(x) ⊗ b(x) p m • ∇ x f + (F ⊥ + F ) • ∇ p f = 0 ( 31
)
where

F ⊥ = -ω(x, p) (p • b(x)) b(x) ∧ (p ∧ b(x)) |p ∧ b(x)| , F = (ω(x, p) |p ∧ b(x)| + qE(t, x) • b(x))b(x)
and

ω(x, p) = |p ∧ b(x)| 2m div x b - (p • b(x)) m ∂ x b b(x) • p |p ∧ b(x)| .
Observe that F ⊥ + F = qb(x) ⊗ b(x)E + ω(x, p) ⊥ p and thus (31) reduces to [START_REF] Bostan | Transport equations with disparate advection fields[END_REF]. Notice that the forces F ⊥ , F may become singular when the momentum p is parallel to the magnetic field. Nevertheless the forces remain bounded around these singularities.

Indeed, since t ∂ x b b(x) = 0 we can write

∂ x b b(x) • p |p ∧ b(x)| = ∂ x b b(x) • b(x) ∧ (p ∧ b(x)) |p ∧ b(x)|
implying that the frequency ω(x, p) remains bounded (and therefore the forces F ⊥ , F as well)

|ω(x, p)| ≤ |p ∧ b(x)| 2m |div x b| + |(p • b(x))| m |∂ x b b(x)| .
It remains to determine the initial condition for (31). For this we multiply (1) with

η(t)ϕ(x, p) where η ∈ C 1 c (R + ) and ϕ ∈ C 1 c (R 3 × R 3 ) ∩ ker T .
We deduce the weak formulation

- R + η ′ (t) R 3 R 3 f ε (t, x, p)ϕ(x, p) dpdxdt -η(0) R 3 R 3 f ε (0, x, p)ϕ(x, p) dpdx - R + η(t) R 3 R 3 f ε (t, x, p) p m • ∇ x ϕ + qE(t, x) • ∇ p ϕ dpdxdt = 0. ( 32
)
Passing to the limit as ε ց 0 one gets

d dt R 3 R 3 f (t, x, p)ϕ(x, p) dpdx = R 3 R 3 f (t, x, p) p m • ∇ x ϕ + qE(t, x) • ∇ p ϕ dpdx and R 3 R 3 f (0, x, p)ϕ(x, p) dpdx = lim tց0 R 3 R 3 f (t, x, p)ϕ(x, p) dpdx = R 3 R 3 f in (x, p)ϕ(x, p) dpdx implying that R 3 R 3 (f in (x, p) -f (0, x, p))ϕ(x, p) dpdx = 0, ∀ ϕ ∈ ker T .
Since we already know that f (0, •, •) ∈ ker T we deduce by Proposition 2.1 that

f (0, x, p) = f in (x, p) =: g in (x, |p ∧ b(x)|, (p • b(x))), (x, p) ∈ R 3 × R 3 (33)
At this stage let us make some comments about the limit model (31), (33). The particles are advected only along the magnetic field lines and consequently there is no current in the orthogonal directions to the magnetic lines

j(t, x) = q R 3 f (t, x, p) p m dp = q R 3 f (t, x, p) p m dp = q R 3 f (t, x, p) (p • b(x)) m dp b(x).
Notice also that the electric field accelerates the particles only along the magnetic field lines. It can be shown that the constraint T f (t) = 0 is propagated by the Vlasov equation (31), i.e., T f (t) = 0 ∀ t > 0 provided that T f (0) = 0. For checking this let us introduce the characteristics (X, P )(s; t, x, p) associated to (31

) dX ds = b(X(s)) ⊗ b(X(s)) P (s) m , dP ds = (F ⊥ + F )(s, X(s), P (s)), (X, P )(t) = (x, p). (34)
A straightforward computation shows that the quantities X(s), |P (s)∧b(X(s))|, (P (s)• b(X(s))) satisfy the characteristic equations ( 27), ( 28), (29) implying that

(X, |P ∧ b(X)|, (P • b(X)) )(0; t, x, p) = (X, R, Z)(0; t, x, |p ∧ b(x)|, (p • b(x))).
Therefore one gets

f (t, x, p) = f (0, X(0; t, x, p), P (0; t, x, p)) = g in (X(0; t, x, p), |P (0; t, x, p) ∧ b(X(0; t, x, p))|, (P (0; t, x, p) • b(X(0; t, x, p)))) = g in ((X, R, Z)(0; t, x, |p ∧ b(x)|, (p • b(x)))) ∈ ker T .
By the above considerations we know that the problems (26), (31) are equivalent.

Nevertheless, for the numerical point of view it is preferable to consider the problem (26) since its phase space (x, r, z) ∈ (R 3 × R + × R) has only 5 dimensions, whereas (31) is posed in a 6 dimensional phase space. At the first sight the resolution of (26) may require some conditions on the boundary r = 0 at any time t > 0. Actually this is not the case, as emphasized in the following proposition.

Proposition 3.1 We denote by µ = µ(x, p) the magnetic moment

µ(x, p) = |p ∧ b(x)| 2 2mB(x) = r 2 2mB(x)
.

If the magnetic field is divergence free, then the magnetic moment is an invariant for (31), resp. (26). In particular the solution of (26) is given by

g(t, x, r, z) =    g in ((X, R, Z)(0; t, x, r, z)), if (t, x, r, z) ∈ R + × R 3 × R ⋆ + × R g in ( X(0; t, x, z), 0, Z(0; t, x, z)), if (t, x, r, z) ∈ R + × R 3 × {0} × R (35)
where ( X, Z) solve

d X ds = Z(s) m b( X(s)), d Z ds = qE(s, X(s)) • b( X(s)), ( X, Z)(t) = (x, z).
Proof. We perform the computations with respect to the coordinates (x, r, z). We have

∂ t + z m b(x) • ∇ x - zr 2m div x b ∂ r + r 2 2m div x b + qE(t, x) • b(x) ∂ z µ = = - r 2 z 2m 2 B 2 (x) (b(x) • ∇ x B + B(x) div x b) = - r 2 z 2m 2 B 2 (x) div x (Bb) = 0.
The invariance of the magnetic moment ensures that R 2 (s; t, x, r, z) > 0 for any s ≥ 0, r > 0 and by the continuity of the application s → R(s; t, x, r, z) we deduce that R(s; t, x, r, z) > 0 for any s ≥ 0, r > 0. In the case r = 0 the same invariance guarantees that R(s; t, x, r, z) = 0 for any s ≥ 0 and thus the characteristic equations for (X, R, Z) reduce to that of ( X, Z). Therefore the solution of ( 26) is given by the formula (35).

) ⊥ p is |p ∧ b(x)| 2 2m div x b b(x) = µ(x, p)B(x) div x b b(x) = -µ(x, p)(∇ x B • b(x)) b(x).
For further computations it is useful to write the equation (31) in conservative form.

By direct calculus we obtain

div x b(x) ⊗ b(x) p m +div p (F ⊥ +F ) = (p • b(x)) 2 |p ∧ b(x)| 2 ∂ x b b(x) • b(x) ∧ (p ∧ b(x)) m (36)
and therefore (31) is equivalent to

∂ t f +div x f b(x) ⊗ b(x) p m +div p (f (F ⊥ +F )) = f (p • b(x)) 2 |p ∧ b(x)| 2 ∂ x b b(x) • b(x) ∧ (p ∧ b(x)) m . (37)
A direct consequence of the above conservative form (with zero average source term)

is the balance of the total energy. Proposition 3.2 We have the kinetic energy balance

d dt R 3 R 3 f (t, x, p) |p| 2 2m dpdx = R 3 E(t, x) • R 3 q f (t, x, p) p m dp dx, t ∈ R + .
In particular si E = -∇ x φ satisfies the Poisson equation -∆ x φ = q ε 0 R 3 f dp then the total energy (kinetic and electrostatic) is conserved

d dt R 3 R 3 f (t, x, p) |p| 2 2m + qφ(t, x) 2 dpdx = 0.
Proof. Notice that the functions f, (p

• b(x)) 2 , |p ∧ b(x)| 2 , |p| 2 belong to ker T and p ∧ b(x) = 0 implying that R 3 R 3 f (t, x, p) |p| 2 2m (p • b(x)) 2 |p ∧ b(x)| 2 ∂ x b b(x) • b(x) ∧ (p ∧ b(x)) m dpdx = 0.
Therefore multiplying (37) by |p| 2 2m yields after integration over

R 3 × R 3 d dt R 3 R 3 f (t, x, p) |p| 2 2m dpdx = R 3 R 3 f (t, x, p) p m • (F ⊥ + F ) dpdx = R 3 E(t, x) • R 3 qf (t, x, p) b(x) ⊗ b(x) p m dp dx = R 3 E(t, x) • R 3 qf (t, x, p) p m dp dx.
The conservation of the total energy follows in standard manner by using the continuity and Poisson equations.

Remark 3.3 The Vlasov equation (31) can be written in conservative form. Indeed, by (36) we have

div x b(x) ⊗ b(x) p m + div p (F ⊥ + F ) = T λ, λ = - (p • b(x)) 2 |p ∧ b(x)| 2 ∂ x b b • p ∧ b(x) mω c
and since T f = 0, the equation (31) is equivalent to

∂ t f + div x f b(x) ⊗ b(x) p m + div p (f (F ⊥ + F -λ ω c (x) p ∧ b(x) )) = 0.
In Proposition 3.1 it was shown that the magnetic moment µ is an invariant for (26). As usual this allows us to reduce (26) to a transport problem depending on one parameter, posed in a 4 dimensional phase space. Indeed, the change of variable

g(t, x, r, z) = k(t, x, µ, z), µ = r 2 2mB(x)
leads to the problem

∂ t k + z m b(x) • ∇ x k + (qE(t, x) -µ∇ x B) • b(x) ∂ z k = 0, (t, x, z) ∈ R + × R 3 × R, µ ∈ R + .
Motivated by such reductions, it is worth searching for other invariants. By direct computation we check the invariance of the energy function.

Asymptotic behaviour

Our goal thereafter will be to give some details about the convergence as ε ց 0 of the solutions (f ε ) ε>0 for (1), ( 2) towards the solution of (31), (33). First we focus on weak convergence results. Secondly we derive strong convergence results for well prepared initial data.

Theorem 4.1 Assume that (f in ε ) ε>0 converges weakly in L 2 (R 3 ×R 3 ) as ε ց 0 to some function f in ∈ L 2 (R 3 × R 3
) and let us denote by (f ε ) ε>0 the sequence of weak solutions of (1) with the initial conditions (f in ε ) ε>0 . We suppose that

E ∈ L ∞ loc (R + ; L ∞ (R 3 )) 3 , b ∈ W 1,∞ (R 3 ) 3 and div x (Bb) = 0. Then (f ε ) ε>0 converges weakly ⋆ in L ∞ (R + ; L 2 (R 3 ×R 3 ))
as ε ց 0 to the weak solution of (31), (33).

Proof. It is a straightforward consequence of the computations in Lemmas 3.2, 3.3.

We only sketch the main steps. Since the characteristic flow associated to the transport

operator p m • ∇ x + q(E(t, x) + p m ∧ B ε ) • ∇ p is measure preserving we deduce that R 3 R 3 |f ε (t, x, p)| 2 dpdx = R 3 R 3 |f in ε (x, p)| 2 dpdx, t ∈ R + , ε > 0
and therefore, after extraction eventually, we can assume that (f ε ) ε>0 converges weakly

⋆ in L ∞ (R + ; L 2 (R 3 × R 3
)) to some function f . Multiplying by ε the weak formulation of (1) written with the test function

η(t)ϕ(x, p), η ∈ C 1 c (R + ), ϕ ∈ C 1 c (R 3 × R 3
) and passing to the limit for ε ց 0 imply that f (t) ∈ ker T , t ∈ R + . Therefore there is a function g = g(t, x, r, z) such that

f (t, x, p) = g(t, x, |p ∧ b(x)|, (p • b(x))), (t, x, p) ∈ R + × R 3 × R 3 . ( 38
)
Choosing now smooth test functions ϕ which belong to the kernel of T yields by letting

ε ց 0 (see (32)) - R + η ′ (t) R 3 R 3 f (t, x, p)ϕ(x, p) dpdxdt -η(0) R 3 R 3 f in (x, p)ϕ(x, p) dpdx - R + η(t) R 3 R 3 f (t, x, p) p m • ∇ x ϕ + qE(t, x) • ∇ p ϕ dpdxdt = 0. ( 39
)
We are done if we prove that the formulation (39) is equivalent to the transport problem (31), (33). For doing this we transform all the integrals in (39) by using the symmetries of f cf. (38) and the properties of the average operator. Since ϕ belongs to the kernel of T , there is a function

ψ such that ϕ(x, p) = ψ(x, |p∧b(x)|, (p•b(x))), (x, p) ∈ R 3 ×R 3 .
It is easily seen, by using cylindrical coordinates along the magnetic axis, that

- R + η ′ (t) R 3 R 3 f ϕ dpdxdt = - R + R 3 R + R g(t, x, r, z)∂ t (ηψ2πr) dzdrdxdt (40)
and -η(0)

R 3 R 3 f in ϕ dpdx = -η(0) R 3 R 3 f in ϕ dpdx = - R 3 R + R g in (x, r, z)(η(0)ψ2πr) dzdrdx. ( 41
)
For the last integral in (39) we appeal to Lemmas 3.2, 3.3.

-

R + η(t) R 3 R 3 f (t, x, p) p m • ∇ x ϕ + qE(t, x) • ∇ p ϕ dpdxdt = - R + η(t) R 3 R 3 f (t, x, p) p m • ∇ x ϕ + qE(t, x) • ∇ p ϕ dpdxdt = - R + η(t) R 3 R 3 f b(x) ⊗ b(x) p m • ∇ x ψ - (p • b(x))|p ∧ b(x)| 2m div x b ∂ r ψ + |p ∧ b(x)| 2 2m div x b ∂ z ψ + q(E(t, x) • b(x))∂ z ψ dpdxdt = - R + R 3 R + R g(t, x, r, z) div x ηψ2πr z m b(x) -∂ r ηψ2πr zr 2m div x b +∂ z ηψ2πr r 2 2m div x b + qE(t, x) • b(x) dzdrdxdt. ( 42
)
Gathering together (40), (41), (42) in (39) yields exactly the weak formulation of (26), (33) and therefore f solves (31), (33). By the uniqueness of the solution of (31), (33) we deduce that all the family (f ε ) ε>0 converges weakly

⋆ in L ∞ (R + ; L 2 (R 3 × R 3 )) towards f .
We inquire now about the strong convergence of the family (f ε ) ε>0 when the initial conditions are well prepared lim εց0

f in ε = f in ∈ ker T strongly in L 2 (R 3 × R 3 ).
These results are definitely much difficult to establish and require smoothness for the solution of the limit model. The key point is the characterization for the solvability of

T u = v (see Proposition 2.2): v ∈ Range T iff v = 0.
Theorem 4.2 Assume that (f in ε ) ε>0 are smooth and converge strongly in

L 2 (R 3 × R 3 ) as ε ց 0 towards some function f in ∈ ker T ∩ C 2 c (R 3 × R 3 ) verifying f in (x, p) = 0, (x, p) ∈ R 3 × R 3 , |p ∧ b(x)| ≤ r in for some r in > 0. We suppose that E ∈ L ∞ loc (R + ; W 1,∞ (R 3 )) 3 , ∂ t E ∈ L ∞ loc (R + ; L ∞ (R 3 )) 3 , b ∈ W 2,∞ (R 3 ) 3 , B ∈ W 1,∞ (R 3 ), div x (Bb) = 0 and inf x∈R 3 B(x) > 0. Then (f ε ) ε>0 con- verges strongly in L ∞ loc (R + ; L 2 (R 3 ×R 3 )) as ε ց 0 to f . In particular if ((f in ε -f in )/ε) ε>0 is bounded in L 2 (R 3 × R 3 ), then ((f ε -f )/ε) ε>0 is bounded in L ∞ loc (R + ; L 2 (R 3 × R 3 )).
Proof. Notice that the characteristics (X, P ) in (34) satisfies

|X(t) -X(0)| ≤ t 0 |P (s)| m ds | |P (t)| -|P (0)| | ≤ |q| t 0 E(s) L ∞ (R 3 ) ds
and therefore it is easily seen that at any time t ∈ R + the solution f of (31), (33) has compact support. Indeed, if

supp f in ⊂ {(x, p) ∈ R 3 × R 3 : |x| ≤ L x , |p| ≤ L p } then for any t ∈ [0, T ], T ∈ R + we have supp f (t, •, •) ⊂ {(x, p) ∈ R 3 × R 3 : |x| ≤ L T x , |p| ≤ L T p }
where

L T p = L p + |q| T 0 E(s) L ∞ (R 3 ) ds, L T x = L x + T m L T p .
Moreover, the regularity of the electro-magnetic field guarantees that for any

T ∈ R + ∇ 2 t,x,p (X, P )(s; t, •, •) ∈ L ∞ (R 3 × R 3 )
uniformly with respect to s, t ∈ [0, T ] and therefore the strong solution f of (31), (33), which is given by

f (t, x, p) = f in ((X, P )(0; t, x, p)) belongs to W 2,∞ ([0, T ] × R 3 × R 3 ).
In particular, the compactness of the support yields

∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f ∈ L ∞ ([0, T ]; W 1,1 ∩ W 1,∞ (R 3 × R 3 )) ∂ t ∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f ∈ L ∞ ([0, T ]; L 2 (R 3 × R 3 )).
The solution f of the limit model ( 31), (33) satisfies

f (0) = f in , ∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f = 0, T f (t) = 0, t ∈ R + .
By Proposition 2.2 there is a unique function h such that

h(t) ∈ D(T ), ∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f + T h(t) = 0, h(t) = 0, t ∈ R + . ( 43
h ∈ L ∞ ([0, T ]; W 1,1 ∩ W 1,∞ (R 3 × R 3 )), ∂ t h ∈ L ∞ ([0, T ]; L 2 (R 3 × R 3 )).
We deduce that

∂ t h + p m • ∇ x h + qE(t, x) • ∇ p h ∈ L ∞ ([0, T ]; L 2 (R 3 × R 3 ))
and thus combining (1), (43) and T f = 0 yields

∂ t (f ε -f -εh) + p m • ∇ x (f ε -f -εh) + qE(t) • ∇ p (f ε -f -εh) + 1 ε T (f ε -f -εh) = -ε ∂ t h + p m • ∇ x h + qE(t) • ∇ p h .
Multiplying by f ε -f -εh and integrating over R 3 × R 3 one gets

1 2 d dt f ε -f -εh 2 ≤ ε ∂ t h + p m • ∇ x h + qE(t) • ∇ p h f ε -f -εh implying that (f ε -f -εh)(t) ≤ f in ε -f in -εh(0) +ε t 0 ∂ t h + p m • ∇ x h + qE(s) • ∇ p h ds, t ∈ R + .
Therefore for any T > 0 there is a constant C T not depending on ε > 0 such that

f ε (t) -f (t) ≤ f in ε -f in + C T ε, t ∈ [0, T ], ε > 0
and our conclusions follow immediately.

Second order gyrokinetic Vlasov equation

In the previous section it was shown that up to O(ε) terms, the particle distributions (f ε ) ε>0 solving (1), ( 2) behave like the solution of the limit model (31) (or equivalently ( 26)) (33).

The main motivation of the guiding-center approximation relies on the confinement properties. We have seen that at the lowest order the particles are advected along the magnetic lines and therefore the plasma remains confined provided that the magnetic field shape is such that the magnetic lines are closed into a bounded domain. Nevertheless, in order to approximate the confinement time we need to compute the drift velocities in the orthogonal directions to the magnetic lines corresponding to the first order correction (X 1 , P 1 ) of the particle dynamics (X ε , P ε ) in [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF]. We inquire about higher order approximation for the Vlasov equation with large magnetic field. Certainly we may expect that the second order approximation follows by similar arguments. Actually this analysis requires much more computations. In this case is convenient to establish first some general results on abstract average operators. It mainly concerns their commutation properties with respect to first order differential operators. Next we will appeal to these results and we will obtain the second order approximation for the Vlasov equation (1).

Average operator revisited

In this section the notation b 0 stands for a given field b

0 : R m → R m satisfying b 0 ∈ W 1,∞ loc (R m ) (44) div y b 0 = 0 ( 45
)
and the growth condition

∃ C > 0 : |b 0 (y)| ≤ C(1 + |y|), y ∈ R m . ( 46
)
Under the above hypotheses the characteristic flow Y = Y (s; y) is well defined

dY ds = b 0 (Y (s; y)), (s, y) ∈ R × R m (47) Y (0; y) = y, y ∈ R m , (48)
and has the regularity Y ∈ W 1,∞ loc (R × R m ). By (45) we deduce that for any s ∈ R, the map y → Y (s; y) is measure preserving

R m θ(Y (s; y)) dy = R m θ(y) dy, ∀ θ ∈ L 1 (R m ).
We have the following standard result concerning the kernel of u → T u = div y (b 0 (y)u(y)).

Proposition 5.1 Let u ∈ L 1 loc (R m ). Then div y (b 0 (y)u(y)) = 0 in D ′ (R m ) iff for any s ∈ R we have u(Y (s; y)) = u(y) for a.a. y ∈ R m .
We denote by T the linear operator defined by T u = div y (b 0 (y)u(y)) for any u in the domain

D(T ) = {u ∈ L 2 (R m ) : div y (b 0 (y)u(y)) ∈ L 2 (R m )}.
Thanks to Proposition 5.1 we have

ker T = {u ∈ L 2 (R m ) : u(Y (s; y)) = u(y), s ∈ R, a.e. y ∈ R m }.
Following the ideas in [START_REF] Bostan | Transport equations with disparate advection fields[END_REF] we introduce the average operator along the measure preserving flow Y .

u = (u -u ) + u , R m (u -u ) u dy = 0.
We will assume that the range of T is closed. This is the case for strongly magnetized plasmas, cf. Proposition 2.2 i.e., Range T = ker • .

(49)

We are looking now for first order differential operators commuting with the average operator. Consider a smooth field c with bounded divergence

c ∈ W 1,∞ loc (R m ), div y c ∈ L ∞ (R m )
and let us denote by Z the flow associated to c (we assume that Z is well defined for any (s, y) ∈ R × R m ). We consider the operator c • ∇ y

D(c•∇ y ) = {u ∈ L 2 (R m ) : div y (cu) ∈ L 2 (R m )}, c•∇ y u = div y (cu)-(div y c)u, u ∈ D(c•∇ y )
where the divergence is understood in distribution sense i.e., there is a function v ∈

L 2 (R m ) such that R m vϕ dy + R m u(c • ∇ y ϕ) dy = 0 for any function ϕ ∈ C 1 c (R m ) or equivalently sup ϕ =0,ϕ∈C 1 c (R m ) R m u(c • ∇ y ϕ) dy ϕ L 2 (R m ) < +∞.
We want to determine the fields c such that c • ∇ y is commuting with • i.e., for any u ∈ D(c • ∇ y ) the average u belongs to D(c • ∇ y ) and

c • ∇ y u = c • ∇ y u .
Notice that the differential operator associated to such a field c leaves invariant the

kernel of T ∀ u ∈ D(c • ∇ y ) ∩ ker T, c • ∇ y u ∈ ker T. (50)
Indeed, for any u ∈ D(c • ∇ y ) ∩ ker T we have u = u and

c • ∇ y u = c • ∇ y u = c • ∇ y u ∈ ker T.
Actually the condition (50) can be written in a simpler form. For any function u ∈ L 2 (R m ) the notation u h stands for the translation u(Z(h; •)). Since div y c is bounded, notice that u h ∈ L 2 (R m ) for any h ∈ R. We appeal to the standard result (see [START_REF] Brezis | Analyse fonctionnelle. Théorie et applications[END_REF],

Proposition IX.3, pp. 153 for similar results).

Lemma 5.1 Let u be a function in L 2 (R m ). Then the following statements are equiv-

alent a) u ∈ D(c • ∇ y ). b) (h -1 (u h -u)) h is bounded in L 2 (R m ).
Moreover, for any u ∈ D(c • ∇ y ) we have the convergence

lim h→0 u h -u h = c • ∇ y u, strongly in L 2 (R m ).
We have the following formula of integration by parts.

Lemma 5.2 For any function u, ϕ ∈ D(c • ∇ y ) we have R m (c • ∇ y ϕ)u dy + R m (c • ∇ y u)ϕ dy + R m (div y c)uϕ dy = 0.
Based on the characterization in Lemma 5.1 we prove that (50) is equivalent to the invariance of ker T by the translations u → u h .

Proposition 5.3 A smooth field c with bounded divergence satisfies (50) iff the translations parallel to c leave invariant the kernel of T i.e.,

∀ h ∈ R, ∀ u ∈ ker T then u h ∈ ker T. (51)
Proof. Assume that (51) holds true and let us consider u ∈ D(c • ∇ y ) ∩ ker T . By Lemma 5.1 we know that

c • ∇ y u = lim h→0 u h -u h strongly in L 2 (R m ).
But for any h = 0 we have (u h -u)/h ∈ ker T and since ker T is closed (because T is closed) we deduce that c • ∇ y u ∈ ker T .

Assume now that (50) holds true and let us establish (51). For the sake of simplicity we assume that the field b 0 possesses a complete family of smooth independent prime

integrals denoted ψ 1 , ..., ψ m-1 i.e., b 0 • ∇ y ψ i = 0, i ∈ {1, 2, ..., m -1}, rank ∂ y ψ(y) = m -1, y ∈ R m
where ψ = t (ψ 1 , ..., ψ m-1 ). Actually this is enough for our purpose since the above hypotheses hold true for strongly magnetized plasmas. Notice that it is sufficient to

consider u = ψ i , i ∈ {1, 2, ..., m -1}. We can write d dh ψ i (Z(h; y)) = (c • ∇ y ψ i )(Z(h; y)). But c • ∇ y ψ i ∈ ker T and therefore c • ∇ y ψ i = v i (ψ 1 (y), ..., ψ m-1 (y)), y ∈ R m for some smooth function v i implying that d dh ψ i (Z(h; y)) = v i (ψ 1 (Z(h; y)), ..., ψ m-1 (Z(h; y))), y ∈ R m .
Similarly we have

d dh ψ i (Z(h; Y (s; y))) = v i (ψ 1 (Z(h; Y (s; y))), ..., ψ m-1 (Z(h; Y (s; y)))), y ∈ R m .
The functions h → ψ(Z(h; y)) and h → ψ(Z(h; Y (s; y))) satisfy the same system of differential equations and the same initial conditions

ψ(Z(0; y)) = ψ(y) = ψ(Y (s; y)) = ψ(Z(0; Y (s; y))).
By the uniqueness of the solution we deduce that

ψ(Z(h; Y (s; y))) = ψ(Z(h; y)), h, s ∈ R
saying that ψ h is constant along the flow Y .

Remark 5.1 In the sequel we will need to pick test functions ϕ ∈ D(c • ∇ y ) ∩ ker T .

When the field b 0 possesses a complete family of smooth independent prime integrals

ψ 1 , ..., ψ m-1 verifying lim |y|→+∞ (|ψ 1 (y)| + ... + |ψ m-1 (y)|) = +∞ (52)
it is easily seen that for any function

v ∈ C 1 c (R m-1 ) the function y → v(ψ 1 (y), ..., ψ m-1 (y)) belongs to C 1 c (R m ) ∩ ker T which is contained in D(c • ∇ y ) ∩ ker T .
Adapting the arguments in the proof of Lemma 5.

D(c • ∇ y ) iff sup ϕ =0,ϕ∈C 1 c (R m )∩ker T R m (c • ∇ y ϕ)u dy ϕ L 2 (R m ) < +∞. (53)
Proof. We only indicate the main lines. The key point is that for functions u ∈ ker T the uniform bound in (53) (when ϕ spans only the constant functions along the flow

Y ) ensures the boundedness of u h -u h h in L 2 (R m
), which implies, by Lemma 5.1, that u ∈ D(c • ∇ y ). Indeed, since by Proposition 5.3 we know that the translation parallel to c leaves invariant ker T , when estimating the L 2 (R m ) norm of u h -u h it is sufficient to integrate against smooth test functions ϕ ∈ C 1 c (R m )∩ker T . After standard computations we obtain a uniform bound for the L 2 norms of u h -u h h provided that (53) holds true.

We are ready now to establish the following result for differential operators leaving invariant the kernel of T .

c • ∇ y u -c • ∇ y u = (u -u ) div y c .
Proof. Let us consider u ∈ D(c • ∇ y ). Using the variational charaterization of the average operator and the integration by parts formula in Lemma 5.2 we can write for

any test function ϕ ∈ C 1 c (R m ) ∩ ker T R m c • ∇ y u ϕ dy = R m (c • ∇ y u) ϕ dy = - R m uϕ div y c dy - R m (c • ∇ y ϕ) u dy. (54)
Since c • ∇ y leaves invariant the kernel of T , the function c • ∇ y ϕ belongs to ker T . We have

R m (c • ∇ y ϕ)u dy = R m (c • ∇ y ϕ) u dy and therefore R m (c • ∇ y ϕ) u dy = - R m (u div y c + c • ∇ y u)ϕ dy implying that sup ϕ =0,ϕ∈C 1 c (R m )∩ker T R m (c • ∇ y ϕ) u dy ϕ L 2 (R m ) ≤ div y c L ∞ (R m ) u L 2 (R m ) + c • ∇ y u L 2 (R m ) .
By Proposition 5.4 we deduce that u ∈ D(c • ∇ y ) and coming back in (54) yields

R m c • ∇ y u ϕ dy = - R m uϕ div y c dy - R m (c • ∇ y ϕ) u dy = R m (c • ∇ y u )ϕ dy + R m ( u -u)ϕ div y c dy = R m (c • ∇ y u )ϕ dy + R m ( u -u)div y c ϕ dy.
Finally one gets for any ϕ ∈ ker T

R m {c • ∇ y u -c • ∇ y u -(u -u )div
(c • ∇ y ϕ )u dy ≤ C ϕ L 2 (R m ) ≤ C ϕ L 2 (R m ) (55)
b i • ∇ y ψ k = δ ik , k ∈ {1, ..., m -1}. ( 57
)
Notice that the differential operators b i • ∇ y leave invariant the kernel of T because for any, let say smooth, u(y) = v(ψ 1 (y), ..., ψ m-1 (y)) ∈ ker T we have

b i • ∇ y u = m-1 k=1 ∂v ∂ψ k b i • ∇ y ψ k = ∂v ∂ψ i (ψ(y)) ∈ ker T. (58)
Therefore, by Proposition 5.5 we have for any i ∈ {1, ..., m -1} and any

u ∈ D(b i • ∇ y ) u ∈ D(b i • ∇ y ), b i • ∇ y u -b i • ∇ y u = (u -u )div y b i .
Obviously the previous statement also holds true for the operator T = b 0 •∇ y . Actually, for any i ∈ {1, ..., m -1}, it is possible to construct a first order differential operator bi • ∇ y commuting with the average operator. Indeed, consider bi = b i -λ i b 0 , where We compute now the fields b i , i ∈ {1, ..., m -1} in the specific case of strongly magnetized plasmas. We have m = 6, y = (x, p), b 0 (x, p) = (0, 0, 0, ω c (x)p ∧ b(x)),

T λ i = div y b i -div y b i , λ i =
T = ω c (x) (p ∧ b(x)) • ∇ p ψ 1 (x, p) = x 1 , ψ 2 (x, p) = x 2 , ψ 3 (x, p) = x 3 , ψ 4 (x, p) = |p ∧ b(x)|, ψ 5 (x, p) = p • b(x) and u (x, p) = 1 2π S(x) u(x, |p ∧ b(x)|ω + (p • b(x))b(x)) dω ( 59
)
where S(x) = {ω ∈ S 2 : b(x) • ω = 0}. Under the assumption inf x∈R 3 B(x) > 0 we know by Proposition 2.2 that (49) holds true. Notice also that (52) is satisfied. We obtain the formula ∇ x,p ψ i = (e i , 0, 0, 0), i ∈ {1, 2, 3}

∇ x,p ψ 4 = - (p • b(x)) |p ∧ b(x)| t ∂ x b p, p -(p • b(x))b |p ∧ b(x)| , ∇ x,p ψ 5 = ( t ∂ x b p, b) b i = e i , - ( t ∂ x b p) i |p ∧ b(x)| ⊥ p , i ∈ {1, 2, 3}, b 4 = 0, 0, 0, p -(p • b(x))b |p ∧ b(x)| , b 5 = (0, 0, 0, b) (60)
where

⊥ p = |p ∧ b(x)| b -(p • b(x)) p-(p•b(x))b |p∧b(x)| . Notice that div p ⊥ p = - (p • b(x)) |p ∧ b(x)| , div x,p b i = ( t ∂ x b p) i |p ∧ b(x)| 2 (p • b(x)), i ∈ {1, 2, 3} div x,p b 4 = 1 |p ∧ b(x)|
, div x,p b 5 = 0.

Multi-scale analysis

The Vlasov equation describing the dynamics of strongly magnetized plasmas (1) is of the form

∂ t u ε + a(t, y) • ∇ y u ε + 1 ε b 0 • ∇ y u ε = 0, t ∈ R + , y ∈ R m (61)
with m = 6, y = (x, p), a(t, x, p) = (p/m, qE(t, x)), b 0 (x, p) = (0, 0, 0, ω c (x)p ∧ b(x)).

Notice that div x,p a = div x,p b 0 = 0. The ansatz u ε = u + εu 1 + ε 2 u 2 + ... leads to the sequence of equations

T u = b 0 • ∇ y u = 0 (62) ∂ t u + a(t, y) • ∇ y u + T u 1 = 0 (63) ∂ t u 1 + a(t, y) • ∇ y u 1 + T u 2 = 0 (64)
. . .

The time evolution equation for the leading order term u comes by applying the average operator to (63) and taking into account that T u 1 ∈ ker • . We obtain Taking into account that u(t) ∈ ker T we have by (58) that b i • ∇ y u(t) ∈ ker T , i ∈ {1, ..., m -1} and b 0 • ∇ y u(t) = 0. Therefore we can write

∂ t u + a(t) • ∇ y u(t) = 0. ( 65
a(t) • ∇ y u(t) = m-1 i=0 α i (t) b i • ∇ y u(t) = m-1 i=1 α i (t) b i • ∇ y u(t). (66)
Actually the formula (66) holds true for any function in the kernel of T saying that any transport operator a • ∇ y reduces, by averaging along the flow of b 0 to another transport operator, denoted a • ∇ y . If a • ∇ y = m-1 i=0 α i b i • ∇ y is a linear combination of differential operators leaving invariant the kernel of T then the effective operator by averaging is

A • ∇ y = a • ∇ y = m-1 i=1 α i b i • ∇ y (67)
and therefore the model for the dominant term u is given by

∂ t u + m-1 i=1 α i b i • ∇ y u = 0.
Remark 5.3 Assume that a, A are smooth and verify

div y a = 0, a • ∇ y u = A • ∇ y u, u ∈ C 1 c (R m ) ∩ ker T.
Integrating with respect to y ∈ R m we deduce that for any

u ∈ C 1 c (R m ) ∩ ker T R m A • ∇ y u dy = R m a • ∇ y u dy = R m a • ∇ y u dy = - R m
u div y a dy = 0.

Therefore we obtain

R m u div y A dy = R m u div y A dy = - R m A • ∇ y u dy = 0
implying that div y A = 0. In particular taking λ 0 = λ 0 (y) such that λ 0 = 0, T λ 0 = div y A (which is possible because div y A = 0) we can replace the averaged transport operator A • ∇ y by an equivalent one

A ′ • ∇ y = A • ∇ y -λ 0 b 0 • ∇ y a • ∇ y u = A • ∇ y u = A ′ • ∇ y u, u ∈ C 1 c (R m ) ∩ ker T such that div y A ′ = div y A -div y (λ 0 b 0 ) = div y A -T λ 0 = 0.
When considering the Vlasov equation (1) and the fields b i , i ∈ {1, 2, 3, 4, 5} in (60) we obtain

α 0 = -m (v ∧ • p) |p ∧ b(x)| 2 , (α 1 , α 2 , α 3 ) = p m , v ∧ = E ∧ b B ( 68
) α 4 = - (p • b(x)) |p ∧ b(x)| ∂ x b : p ⊗ p m + qE • p -(p • b(x))b |p ∧ b(x)| , α 5 = ∂ x b : p ⊗ p m + qE • b. (69)
The average of the above coefficients are

α 0 = 0, (α 1 , α 2 , α 3 ) = (p • b(x)) m b (70)
α 4 = - (p • b(x))|p ∧ b(x)| 2m div x b, α 5 = |p ∧ b(x)| 2 2m div x b + qE • b (71)
and therefore the average of the transport operator p

m • ∇ x + qE • ∇ p is p m • ∇ x + qE • ∇ p = (p • b(x)) m b • ∇ x + (q(E • b)b + ω(x, p) ⊥ p) • ∇ p with ω(x, p) = |p ∧ b(x)| 2m div x b - (p • b(x)) m ∂ x b b • p |p ∧ b(x)| .
We have retrieved the limit model (31).

((p • b) b/m, q(E • b) b + ω(x, p) ⊥ p) if E(t, x) = -∇ x φ(t, x).
In the sequel we intend to take into account the first order corrections u 1 , that is, write a model whose solution coincides with u + εu 1 up to O(ε 2 ). By using the orthogonal

decomposition in L 2 (R m ) u 1 = v 1 + w 1 , T v 1 = 0, w 1 = 0
we obtain from (63), (65)

-T w 1 = ∂ t u + a • ∇ y u = ∂ t u + a • ∇ y u -∂ t u + a • ∇ y u = m-1 i=1 (α i -α i )b i • ∇ y u.
Therefore the zero average contribution of u 1 is given by

w 1 = -Ã • ∇ y u, Ã = m-1 i=1 β i b i (72)
where for any i ∈ {1, 2, ..., m -1}, β i solves

T β i = α i -α i ∈ ker • , β i = 0.
Applying the average operator in (64) yields

∂ t v 1 + A • ∇ y v 1 + a • ∇ y w 1 = 0 (73)
and combining (65), (73) leads to

∂ t (u + εu 1 ) + A • ∇ y (u + εu 1 ) + ε( a • ∇ y w 1 -∂ t w 1 -A • ∇ y w 1 ) = 0. ( 74
)
Replacing the time derivative in terms of space derivatives thanks to ∂ t u = -A • ∇ y u we transform the time derivative of w 1 as follows

-∂ t w 1 = ∂ t à • ∇ y u -à • ∇ y (A • ∇ y u).
The equation (74) can be written in the form

∂ t (u + εu 1 ) + A • ∇ y (u + εu 1 ) + ε(∂ t à • ∇ y u + R(u)) = 0 (75)
where

R(u) = a • ∇ y w 1 -A • ∇ y w 1 -Ã • ∇ y (A • ∇ y u).
We need to express R(u) in terms of u in view of (72). Eventually R(u) contains second order derivatives of u. We will see that all second order derivatives cancel. Moreover we will show that R(u) reduces to a first order differential operator. Let us split the computations into several steps.

Proposition 5.6 Assume that div y a = 0. For any function u ∈ ker T we have

a • ∇ y w 1 = m-1 j=1 m-1 i=1 {b i • ∇ y β i T β j + div y b i β i T β j } b j • ∇ y u.
Proof. By Proposition 5.5 we can write

a • ∇ y w 1 = m-1 i=0 α i b i • ∇ y w 1 = m-1 i=0 b i • ∇ y (α i w 1 ) - m-1 i=0 (b i • ∇ y α i )w 1 = m-1 i=0 b i • ∇ y α i w 1 - m-1 i=0 (α i w 1 -α i w 1 )div y b i - m-1 i=0 (b i • ∇ y α i )w 1 = m-1 i=0 b i • ∇ y α i w 1 + m-1 i=0 α i w 1 div y b i -w 1 div y a = m-1 i=1 b i • ∇ y α i w 1 + m-1 i=1 α i w 1 div y b i .
In the last equality we have taken into account that α 0 w 1 ∈ ker T and div y b 0 . Since w 1 = 0 we have for any i ∈ {1, ..., m -1}

α i w 1 = w 1 T β i = -β i T w 1 = m-1 j=1 β i T β j b j • ∇ y u.
Notice that the matrix ( β i T β j ) i,j is skew-symmetric

β i T β j + β j T β i = T (β i β j ) = 0, i, j ∈ {1, ..., m -1}.
Another key point is that for any i, j ∈ {1, ..., m -1} the operators b i • ∇ y , b j • ∇ y are commuting on functions of ker T . Indeed if u = v(ψ 1 , ..., ψ m-1 ) ∈ ker T we have cf.

(

) b i • ∇ y u = ∂ ψ i v(ψ(y)), b j • ∇ y u = ∂ ψ j v(ψ(y)). It follows that b i • ∇ y (b j • ∇ y u) -b j • ∇ y (b i • ∇ y u) = ∂ 2 ψ i ψ j v(ψ(y)) -∂ 2 ψ j ψ i v(ψ(y)) = 0. 58
Combining the skew-symmetry of ( β i T β j ) i,j and the symmetry of (b

i •∇ y (b j •∇ y u)) i,j we deduce that m-1 i=1 b i •∇ y α i w 1 = m-1 i=1 b i •∇ y m-1 j=1 β i T β j b j • ∇ y u = m-1 j=1 m-1 i=1 b i • ∇ y β i T β j b j •∇ y u.
Therefore the term a • ∇ y w 1 can be written

a • ∇ y w 1 = m-1 j=1 m-1 i=1 {b i • ∇ y β i T β j + div y b i β i T β j } b j • ∇ y u.
Proposition 5.7 For any function u ∈ ker T we have

-A • ∇ y w 1 -Ã • ∇ y (A • ∇ y u) = [A, Ã] • ∇ y u
where [A, Ã] stands for the Poisson bracket between the fields

A = m-1 i=1 α i b i and à = m-1 i=1 β i b i .
Proof. Obviously we have

-A • ∇ y w 1 -Ã • ∇ y (A • ∇ y u) = A • ∇ y ( Ã • ∇ y u) -Ã • ∇ y (A • ∇ y u) = [A, Ã] • ∇ y u.
Proposition 5.8 Assume that div y a = 0. Then for any function u ∈ ker T we

R(u) = m-1 j=1 m-1 i=1 {b i • ∇ y β i T β j + div y b i β i T β j } b j • ∇ y u + [A, Ã] • ∇ y u. (76) Remark 5.5 If I(t) ∈ ker T is a prime integral for a(t) = m-1 i=0 α i (t)b i therefore I(t) is also a prime integral for A(t) = m-1 i=1 α i (t) b i cf. Remark 5.

and thus the corresponding function w

1 vanishes, because -T w 1 = a(t) • ∇ y I(t) -A(t) • ∇ y I(t) = 0.
We deduce by Proposition 5.8 that I(t) is a prime integral for R i.e., R(I(t)) = 0.

Coming back in (75) we deduce that u + εu 1 solves

∂ t (u + εu 1 ) + A • ∇ y (u + εu 1 ) + ε[∂ t à • ∇ y (u + εu 1 ) + R(u + εu 1 )] = ε 2 [∂ t à • ∇ y u 1 + R(u 1 )].
Therefore we expect that the function ũε solving ∂ t ũε + (A + εA 1 ) • ∇ y ũε = 0 (77)

A 1 = ∂ t à + m-1 j=1 m-1 i=1 {b i • ∇ y β i T β j + div y b i β i T β j } b j + [A, Ã]
will be a second order approximation for u + εu 1 ũε = u + εu 1 + O(ε 2 ).

Consequently, motivated by the ansatz u ε = u + εu 1 + ε 2 u 2 + ..., the solution of (77)

will approximate the solution of (61) up to a second order term ũε = u ε + O(ε 2 ).

Derivation of the gyrokinetic Vlasov equation

This section is devoted to the explicit computation of the second order model (77) which corresponds to the Vlasov equation ( 1) with strong magnetic field. We already know that the transport operator A • ∇ x,p has the form

A • ∇ x,p = (p • b) m b • ∇ x + (q(E • b) b + ω(x, p) ⊥ p) • ∇ p .
It remains to identify the transport operator A 1 ∇ x,p . We need the following lemma, whose proof is left to the reader. for some r in > 0. We denote by f = f (t, x, p) the solution of the limit model

∂ t f + b(x) ⊗ b(x) p m • ∇ x f + qb(x) ⊗ b(x)E + ω(x, p) ⊥ p • ∇ p f = 0, (p ∧ b) • ∇ p f = 0
with the initial condition f (0, x, p) = f in (x, p), (x, p) ∈ R 3 × R 3 . For any t ∈ R + we consider h(t) ∈ D(T ) the unique function satisfying

∂ t f + p m • ∇ x f + qE(t, x) • ∇ p f + T h(t) = 0, h(t) = 0.
Therefore for any T > 0 the function h belongs to W 1,∞ ([0, T ] × R 3 × R 3 ) and has compact support. In particular

∂ t h + p m • ∇ x h + qE(t, x) • ∇ p h ∈ L ∞ ([0, T ]; L 2 (R 3 × R 3 )).
Proof. As in the proof of Theorem 4.2, the regularity of the electro-magnetic field guarantees that f ∈ W 2,∞ ([0, T ] × R 3 × R 3 ). Moreover, the compactness of the support of f in implies the compactness of the support of f (t), uniformly for t ∈ [0, T ]. We have -T h(t) = We need to compute the matrix ( β i T β j ) 1≤i,j≤5 . The following formula will be used. Finally combining Propositions 5.9, 5.12 we obtain Theorem 5.2 The second order model (77) which corresponds to the Vlasov equation

(1) with strong magnetic fields is

∂ t f ε + (A x + εA 1 x ) • ∇ x f ε + (A p + εA 1 p ) • ∇ p f ε = 0 ( 83
)
where

A x = (p • b) m b, A p = q(E • b) b + ω(x, p) ⊥ p, ε Ãx = - p ∧ b mω ε c ε Ãp = m v ε ∧ • p |p ∧ b| p -(p • b)b |p ∧ b| + [ t ∂ x b p • (p ∧ b) + ω c (∂ x b : T -1 (p ⊗ p -p ⊗ p )] ⊥ p mω ε c |p ∧ b| εA 1 x = v ε ∧ + v ε GD + v ε CD + v ε RD + [A, ε Ã] x εA 1 p = m ∂ t v ε ∧ • p |p ∧ b| p -(p • b)b |p ∧ b| - t ∂ x b p |p ∧ b| • (v ε ∧ + v ε GD + v ε CD + v ε RD ) ⊥ p + (p • b)|p ∧ b| 2m 2 ω ε c B ε div x b ∧ ∂ x b b B ε ⊥ p + (∂ x b b • v ε ∧ ) 3 2 (p • b) b - p 2 -(div x v ε ∧ ) p -(p • b)b 2 + [A, ε Ã] p and B ε = B ε , ω ε c = qB ε m v ε ∧ = E ∧ b B ε , v ε GD = |p ∧ b| 2 2m 2 ω ε c b ∧ ∇ x B ε B ε , v ε CD = (p

  • p |p ∧ b(x)| , p ∧ b(x) = 0 and for any (x, p) such that p ∧ b(x) = 0 the symbol ⊥ p stands for the orthogonal momentum to p in the plan determined by b(x) and p such that its coordinate along b(x) is positive ⊥ p = |p ∧ b(x)| b(x) -(p • b(x)) p -(p • b(x))b(x) |p ∧ b(x)| .

  with P (s) and b(x) we deduce that |P (s)| 2 = |p| 2 and b(x) • P (s) = b(x) • p implying also that |b(x) ∧ P (s)| = |b(x) ∧ p|. If the initial conditions satisfy b(x) ∧ p = 0 then clearly P (s; x, p) = (b(x) • p) b(x). If b(x) ∧ p = 0 we consider the positive oriented basis of R 3

  = (p • b(x)) and finally P (s; x, p) = cos(ω c (x)s) b(x) ∧ (p ∧ b(x)) + sin(ω c (x)s) p ∧ b(x) + (b(x) • p) b(x).

Remark 3 . 1 Remark 3 . 2

 3132 If a charged particle situated at the point x has momentum parallel to the magnetic field, i.e., p ∧ b(x) = 0, then at any time s the particle remains to the same magnetic line and its coordinates in the phase space satisfy dX ds = b(X(s)) • P (s) m b(X(s)), d ds (b(X(s))•P (s)) = q E(s, X(s))•b(X(s)), P (s)∧b(X(s)) = 0. We recognize the expression of the diamagnetic force acting in the parallel direction to the magnetic field. The parallel component of the force ω(x, p

Proposition 3 . 3

 33 Assume that the electric potential is stationary i.e., ∂ t φ = 0. Then the energy e = |p| 2 2m + qφ(x) = r 2 +z 2 2m + qφ(x) is an invariant for (31), resp. (26).

) By Remark 2 . 1 ,

 21 Corollary 2.1 and Proposition 5.10 the corrector function h(t, •, •) has compact support in R 3 × R 3 , uniformly with respect to t ∈ [0, T ], and the regularity

Proposition 5 . 2

 52 For any function u ∈ L 2 (R m ) the averages (T -1 T 0 u(Y (s; •)) ds) T >0 , (T -1 0 -T u(Y (s; •)) ds) T >0 converge strongly in L 2 (R m ) as T → +∞ towards some function denoted u ∈ ker T . The average operator u → u is linear continuous on L 2 (R m ). Moreover it coincides with the orthogonal projection on the kernel of T i.e., u ∈ ker T : R m (u -u )ϕ(y) dy = 0, ϕ ∈ ker T. As in Section 2 we have ker • = ker Proj kerT = ⊥ (ker T ) = ⊥ (ker T ⋆ ) = Range T and for any u ∈ L 2 (R m ) we have the orthogonal decomposition

Proposition 5 . 4

 54 1 yields a similar characterization for the elements of D(c • ∇ y ) ∩ ker T . Consider a smooth field c = c(y) with bounded divergence such that c • ∇ y leaves invariant the kernel of T and a function u ∈ ker T . Then u belongs to

Proposition 5 . 5

 55 Consider a smooth field c = c(y) with bounded divergence such that the operator c • ∇ y leaves invariant the kernel of T . Then for any function u ∈ D(c • ∇ y ) the average u belongs to D(c • ∇ y ) and

Remark 5 . 2

 52 y c }ϕ dy = 0 and since c • ∇ y u ∈ ker T (because c • ∇ y leaves invariant the kernel of T ) we deduce that c • ∇ y u -c • ∇ y u = (u -u ) div y c . We may expect a simpler justification for Proposition 5.4 by extending the uniform bound in (53) to any smooth function in L 2 (R m ), not necessarily in ker T . Indeed for any ϕ ∈ C 1 c (R m ) we can use the orthogonal decomposition in L 2 (R m ) ϕ = ϕ + (ϕ -ϕ ) then try to write for some constant C R m

(

  since ϕ ∈ ker T ) by expecting that c • ∇ y (ϕ -ϕ ) = 0, which can be motivated by the fact that ϕ -ϕ = 0. (56) Actually (56) is not valid because it is not of all clear that c • ∇ y leaves invariant the kernel of the average operator. Indeed, by Proposition 5.5 we have for any zero average function θ c • ∇ y θ = -θ div y c which clearly says that the last assertion is false at least when div y c is not constant along the flow Y . Conversely, assume that c • ∇ y leaves invariant the kernel of T and that div y c is constant along the flow of b 0 . By Proposition 5.5, for any u ∈ D(c • ∇ y ) the average u belongs to D(c • ∇ y ) and c • ∇ y u -c • ∇ y u = (u -u ) div y c = u -u div y c = 0. For any i ∈ {1, 2, ..., m -1} consider the field b i orthogonal to b 0 and satisfying

  0. It easily seen that bi • ∇ y leaves invariant the kernel of T and that div y bi is constant along the flow of b 0 div y bi = div y b i -T λ i = div y b i . Therefore Theorem 5.1 implies that bi •∇ y and • are commuting for any i ∈ {1, ..., m-1}.

)

  Consider now the m -1 fields b 1 , ..., b m-1 orthogonal to b 0 and satisfying (57). Let us denote by α 0 , α 1 , ..., α m-1 the coordinates of a in the basis b 0 , b 1 , ..., b m-1 a(t, y) = m-1 i=0 α i (t, y)b i (y).

Remark 5 . 4

 54 It is also possible to determine the average transport operator A • ∇ y = a • ∇ y by imposing (66) for any prime integral ψ i , i ∈ {1, 2, ..., m -1} together with the condition A • b 0 = 0. It is easily seen that any prime integral for a which belongs to the kernel of T is also a prime integral for A. Indeed, if a • ∇ y I = 0 then A • ∇ y I = a • ∇ y I = 0. In particular |p| 2 /2m + qφ(t, x) ∈ ker T is a prime integral for

Lemma 5 . 3 Proposition 5 . 9 2 .

 53592 We have the formulaT -1 (p -p ) = -p ∧ b ω c (x)(78)T -1 (p ⊗ p -p ⊗ p ) = -The transport operator à • ∇ x,p = Ãx • ∇ x + Ãp • ∇ p is given by Ãx = -p ∧ b mω c Ãp = m(v ∧ • p) |p ∧ b(x)| p -(p • b(x))b |p ∧ b(x)| +[ t ∂ x b p•(p∧b)+ω c (∂ x b : T -1 (p⊗p-p ⊗ p )] ⊥ p mω c |p ∧ b| .In particular we have∂ t à • ∇ x,p = m(∂ t v ∧ • p) |p ∧ b(x)| p -(p • b(x))b |p ∧ b(x)| • ∇ p , v ∧ = E ∧ b B .(80)Proof. Using (78) and the expressions of the coefficients α i , α i , i ∈ {1, ...5} computed before cf. (68), (69), (70), (71) we obtain(β 1 , β 2 , β 3 ) = T -1 p -p m = -p ∧ b mω c β 4 = -(p • b(x)) |p ∧ b(x)| ∂ x b : T -1 p ⊗ p -p ⊗ p m + m v ∧ • p |p ∧ b(x)| β 5 = ∂ x b : T -1 p ⊗ p -p ⊗ p m .Therefore thanks to (60) we deduce thatÃx = -p ∧ b mω c (x) Ãp = t ∂ x b p • (p ∧ b) + ω c (∂ x b : T -1 (p ⊗ p -p ⊗ p ) ⊥ p mω c |p ∧ b| +m v ∧ • p |p ∧ b| p -(p • b)b |p ∧ b| and(80) follows.Based on the above considerations we complete now the regularity result used in the proof of Theorem 4.Proposition 5.10 We suppose thatE ∈ L ∞ loc (R + ; W 1,∞ (R 3 )) 3 , ∂ t E ∈ L ∞ loc (R + ; L ∞ (R 3 )) 3 , b ∈ W 2,∞ (R 3 ) 3 , B ∈ W 1,∞ (R 3 ), div x (Bb) = 0 and inf x∈R 3 B(x) > 0. Let f in ∈ ker T ∩ C 2 c (R 3 × R 3 ) verifying f in (x, p) = 0, (x, p) ∈ R 3 × R 3 , |p ∧ b(x)| ≤ r in

5 i=1( 5 i=1T 5 i=1β

 555 α i -α i ) b i • ∇ x,p f (t) = β i b i • ∇ x,p f (t) saying that -h(t) = i b i • ∇ x,p f (t) = à • ∇ x,p f (t).Obviously h has compact support and by Lemma 5.3, Proposition 5.9 we deduce that the components of the field à are locally Lipschitz in{(t, x, p) ∈ [0, T ] × R 3 × R 3 : |p ∧ b(x)| ≥ r in inf B/ sup B}.Our conclusion follows immediately by observing that the invariance of the magnetic moment µ= |p ∧ b| 2 /2mB ensures supp f (t) ⊂ (t, x, p) ∈ [0, T ] × R 3 × R 3 : |p ∧ b(x)| ≥ r in inf B sup B , t ∈ R + .

Lemma 5 . 4 2 2m 2 2 2m 2 ωc b ∧ ∂ x b b |p∧b| 2 t

 542222 For any k ∈ {1, 2, 3} we havep k p ⊗ p = (p•b)b k (p • b) 2 -3 2 |p ∧ b| 2 b⊗b+ (p • b)|p ∧ b| 2 2 [b k I +b⊗e k +e k ⊗b].(81)Proposition 5.11 The elements of the matrix (β i T β j ) 1≤i,j≤5 are given by ωc M [b] -|p∧b| 2 (v ∧ + v CD ) (p•b)|p∧b| (v ∧ + v CD ) 0 (p•b)|p∧b| 2B E • (b ∧ ∂ x b b) -t (p•b)|p∧b| 2 2m 2 ωc b ∧ ∂ x b b -(p•b)|p∧b| 2B E • (b ∧ ∂ x b b)where M [b] is the matrix of the linear map p → b ∧ p M [b] = and v ∧ = E∧b B , v CD = (p•b) 2 m 2 ωc b ∧ ∂ x b b.

• b) 2 m 2 ω ε c b∧∂ x b b, v ε RD = |p ∧ b| 2 2m 2

 2 ω ε c (b•rot x b)b.

  • . The key point here is the Poincaré inequality Proposition 2.2 We assume that inf x∈R 3 B(x) > 0. Then T restricted to ker • is one to one map onto ker • . Its inverse belongs to L(ker • , ker • ) and we have the By the previous computations we know that Range T ⊂ ker • . Assume now that u ∈ D(T ) ∩ ker • such that T u = 0. Since • = Proj ker T we have u = u = 0 saying that T | ker • is injective. Consider now v ∈ ker • and let us prove that there is

	Poincaré inequality
	u ≤	2π |ω 0 |	T u , ω 0 =	q m	inf x∈R 3 B(x) = 0	(17)
	for any u ∈ D(T ) ∩ ker • .
	Proof.

We formulate now necessary and sufficient conditions for fields to be commuting with the average operator.

Therefore for any ϕ ∈ ker T we have

We deduce that ϕ div y c is constant along the flow of b 0 for any function ϕ ∈ ker T . We are done if we are able to construct functions ϕ ∈ ker T with support including balls in R m with radii as large as we want. For example, consider

Clearly the condition (52) implies that for any R > 0, ϕ R ∈ C 1 c (R m ) ∩ ker T and for any y ∈ B(0, R) we have ϕ R (y) = 0 saying that B(0, R) ⊂ supp ϕ R .

Proof. Recall that the matrix ( β i T β j ) 1≤i,j≤5 is skew-symmetric and β i T β j = β i α j , i, j ∈ {1, ..., 5}.

After computations we deduce that the matrix β ⊗ α is given by (82).

Thanks to Proposition 5.6 we obtain the following explicit expression for the term a • ∇ y w 1 . Proposition 5.12 For any function u ∈ ker T we have

where v ∧ , v GD , v CD , v RD are the electric cross field drift, the magnetic gradient drift, the magnetic curvature drift and the magnetic rotational drift respectively

Proof. The entries of the matrix ( β i T β j ) 1≤i,j≤5 belong to the kernel of T . Thus the terms b i • ∇ x,p β i T β j can be easily computed since b i • ∇ x,p = ∂ ψ i for any i ∈ {1, ..., 5} cf. (58). Notice also that div x,p b i = 0 for i ∈ {1, 2, 3, 5} and div x,p b

After computations one gets

and the conclusion follows.


Sébastien Lengagne
email: lengagne@lirmm.fr

Nacim Ramdani
email: ramdani@lirmm.fr

Philippe Fraisse
email: fraisse@lirmm.fr

A new method for generating safe motions for humanoid robots

Keywords: Motion planning, Semi Infinite Programming, Discretization, Interval Analysis, Constraints

This paper introduces a new method for planning safe motions for complex systems such as humanoid robots. Motion planning consists on finding the best joint trajectories. By using trajectory parameterization, the motion planning problem can be seen as a Semi-Infinite Programming problem (SIP) since it involves a finite number of parameters over an infinite set of constraints. Most methods solve the SIP problem by transforming it into a finite programming one by using a discretization over a prescribed grid. We show that this approach is risky because it can lead to motions which violate one or several constraints. Then we introduce our new method for planning safe motions. It uses Interval Analysis techniques in order to achieve a safe discretization of the constraints. We show how to implement this method and use it with state-of-the-art constrained optimization packages. Then, we illustrate its capabilities for planning safe motions for the HOAP-3 humanoid robot.

I. INTRODUCTION

Motion planning is a classical topic in robotics research. For complex systems such as humanoid robots, motion planning deals with a great number of constraints, which increases the computation time. Therefore, motion planning is usually done off-line to create a database of benchmark motions [START_REF] Denk | Synthesis of a Walking Primitive Database for a Humanoid Robot using Optimal Control Techniques[END_REF]. Then a simple on-line control loop can track one of these motions.

Motion planning includes the problem of digital actors' locomotion [START_REF] Yoshida | Humanoid motion planning for dynamic tasks[END_REF], kick motion generation on HRP-2 robot [START_REF] Miossec | Development of a software for motion optimization of robots -application to the kick motion of the hrp-2 robot[END_REF], computing a manipulator robot's trajectory [START_REF] Piazzi | Global minimum-jerk trajectory planning of robot manipulators[END_REF] or smoothing pre-calculated motions [START_REF] Suleiman | On humanoid motion optimization[END_REF]. In these works, the planned motions minimize a cost function and validate sets of equality and inequality constraints. With a joint trajectory parameterization, we transform the motion planning problem into a Semi-Infinite Programming problem (SIP), since it involves a finite number of parameters over an infinite continuous set of constraints.

To solve a SIP problem, the set of continuous inequality constraints are usually discretized by picking up several values over a given grid. Therefore, the obtained motions will satisfy the constraints only for the grid nodes. However, between two nodes the retained motion may violate some constraints that can have disastrous consequences on the integrity of the systems.

This paper presents a new method for planning safe motions, i.e. motions which ensure that the inequality constraints remain satisfied all over the motion duration. Our method uses the same optimization algorithms as classical one but replace the time-point discretization by a safe discretization that computes the constraints over time-intervals using Interval Analysis [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF]. Interval Analysis has already been used in order to solve in a guaranteed way the problems of self-collision avoidance and prevention for the arms of a 2-degrees of freedom robot [START_REF] Fang | Dynamic interference avoidance of 2-DOF robot arms using interval analysis[END_REF] or to solve the problem of finding collision-free paths [START_REF] Jaulin | path planning using intervals and graphs[END_REF].

A preliminary version of our safe discretization method was tested successfully on a two degrees of freedom pendulum where an optimal one-step motion was generated [START_REF] Lengagne | Guaranteed computation of constraints for safe path planning[END_REF]. In the present paper, we address motion planning issues for a more complex system with six degrees of freedom. To be able to yield results in tractable CPU time, our method has been significantly improved. To illustrate the capabilities of our new method, we generate optimal step motions for the HOAP-3 humanoid robot in the sagital plane while considering two different cost functions and also using two optimization algorithms: C-FSQP [START_REF] Lawrence | User's Guide for CFSQP Version 2.5: A C Code for Solving (Large Scale) Constrained Nonlinear (Minimax) Optimization Problems, Generating Iterates Satisfying All Inequality Constraints[END_REF] and IPOPT [START_REF]Introduction to IPOPT : a tutorial for downloading, installing and using IPOPT[END_REF].

This paper is organized as follows: in section II we introduce the 2-D model of the humanoid robot used to generate motion in the sagital plane. Section III presents how the motion planning problem is transformed into a Semi-Infinite Programming problem, by a joint trajectory parameterization and emphasizes how risky the usual way of constraint discretization is. Section IV introduces the main ideas of our safe motion planning method and its practical implementation via Interval Analysis. In section V, the comparison between the classical and the safe motion planning methods is exposed for the two analyzed cost functions.

II. MODELING

A. 6 dof serial chain

Let us consider the motions of the HOAP-3 humanoid robot in the sagital plane and suppose that the balance in the frontal plane may be controlled using, for instance, the method presented in [START_REF] Fraisse | Towards dynamic balance control of humanoid robots by using com and zmp[END_REF]. As a consequence, we can reduce the 3-D humanoid robot into a 2-D plane system, assuming that the upper parts of the body are equivalent to the chest (Cf. Fig 1). Furthermore, we will consider a serial chain starting from the fixed foot and ending at the flying foot with 6 degrees of freedom (2 ankles, 2 knees, 2 hips).

B. Joint trajectory

We define a motion via the vector X = [τ, p 1 , p 2 , . . . , p 6 ] where τ is the motion duration and p i the coefficients of the weighted sum of B-spline functions which model the i th joint position trajectory q i (t), as follows:

q i (t) = N s ∑ j=0 p i, j × B j (t) (1)
The joint velocity and acceleration are obtained by differentiating 1.

qi (t) = N s ∑ j=0 p i, j × Ḃ j (t) (2)
qi (t) = N s ∑ j=0 p i, j × B j (t) (3)

C. Dynamic modeling

The dynamic model equation allows to compute the joint torques Γ(t) knowing the joint angle, velocity and acceleration vectors. This paper deals with the optimization of motions in the sagital plane, therefore we use the classical Newton-Euler method [START_REF] Khalil | Modeling, Identification & Control of Robots[END_REF] for computing the joint torques of 2-D robots:

Γ(t) = NE(q(t), q(t), q(t),t) (4)

D. Balance

The computation of the Zero Moment Point (ZMP) gives information about the balance of humanoid robots. The ZMP is defined in [START_REF] Vukobratovic | Contribution to the synthesis of biped gait[END_REF] as a point, on the contact surface, where total inertia force is equal to zero. If this point stays within the base of support, the robot will keep its balance. The ZMP depends on the joint angle, velocity and acceleration. ZMP(t) = f (q(t), q(t), q(t),t).

(5)

III. MOTION PLANNING

A. Semi-Infinite Programming

For motion planning, one usually solves a constrained optimization problem in order to find out the optimal joint trajectory q(t) which: minimizes F(q(t)) (6) subject to ∀i, ∀t ∈ [0, T ] g i (q(t)) ≤ 0 (7) and ∀ j h j (q(t)) = 0 [START_REF] Jaulin | path planning using intervals and graphs[END_REF] where F denotes the cost or objective function, g i the set of inequality constraint functions and h i the set of equality constraint functions. This is a classical optimal control problem where there are a continuous function: q(t) and a set of continuous inequality constraints g i (q(t)). Therefore we have to deal with an infinite dimensional problem. The B-spline parameterization allow to have a finite number of value to optimized and transform the problem into finding out the optimal parameter vector X which:

minimizes F(X) (9) subject to ∀i, ∀t ∈ [0, T ] g i (X,t) ≤ 0 (10) and ∀ j h j (X) = 0 ( 11
)
This problem is defined as a Semi infinite Programming problem (SIP) , because there is a finite number of constraints which must be satisfied over a continuous time that can be expressed as an infinite number of instant t [START_REF] Hettich | Semi-infinite programming: theory, methods, and applications[END_REF].

B. Cost function

The choice of the cost function F(X) for motion planning must take into account the features of the robot and the desired application. Some authors minimize motion duration [START_REF] Piazzi | Global minimum-time trajectory planning of mechanical manipulators using interval analysis[END_REF] or jerk [START_REF] Piazzi | Global minimum-jerk trajectory planning of robot manipulators[END_REF] for robot manipulators. Some others minimize the energy consumption taking into account the parameters of the motors (friction, ...) [START_REF] Miossec | Development of a software for motion optimization of robots -application to the kick motion of the hrp-2 robot[END_REF] or a biological inspired cost function: the torque change [START_REF] Uno | Formation and control of optimal trajectory in human multijoint arm movement[END_REF].

In this paper we investigate two different cost functions:

• the motion duration F(X) = T 0 1dt (Cf. Sections III-F and IV-D) • the mechanical energy consumption which is the integral of the Euclidean norm of joint torques

F(X) = T 0 Γ T Γdt (Cf. Section V)

C. Equality constraint functions

The set of the equality constraint functions h j (X) allows to define the motion. These functions usually correspond to constraints on some system state variables at given time instants such as the beginning or the end of a motion. In our case, we consider six equality constraints which are the position of the flying foot at the beginning and at the end of the motion.

D. Inequality constraint functions

The set of the inequality constraints g i (X) is build from the physical limits of the system. These inequalities must be satisfied over the whole motion duration. In these experiments we consider limitations on the joint position, velocity and torque values and also on the ZMP location in order to ensure the robot balance. The set of inequality constraint functions is as follows:

       q ≤ q(t) ≤ q q ≤ q(t) ≤ q Γ ≤ Γ(t) ≤ Γ ZMP ≤ ZMP(t) ≤ ZMP (12)

E. Classical constraint discretization

Classical optimization algorithms use a finite number of discrete constraints. Hence constraints must be discretized. We present the classical way of discretizing them and we emphasize the fact that some constraints can be violated.

In the context of motion planning, discretization usually consists in picking up several time points in a grid for computing the constraint functions [START_REF] Stryk | Numerical solution of optimal control problems by direct collocation[END_REF], [START_REF] Hettich | Semi-infinite programming: theory, methods, and applications[END_REF]. Therefore, the inequality constraints in [START_REF] Lawrence | User's Guide for CFSQP Version 2.5: A C Code for Solving (Large Scale) Constrained Nonlinear (Minimax) Optimization Problems, Generating Iterates Satisfying All Inequality Constraints[END_REF] are replaced by:

∀i, ∀t k ∈ T g i (X,t k ) ≤ 0 ( 13
)
where

T = {t 0 = 0,t 1 , ...,t N-1 ,t N = T }.
Consequently, the continuous problem [START_REF] Lawrence | User's Guide for CFSQP Version 2.5: A C Code for Solving (Large Scale) Constrained Nonlinear (Minimax) Optimization Problems, Generating Iterates Satisfying All Inequality Constraints[END_REF] where appears ∀t ∈ [0, T ] becomes a discrete one: ∀t k ∈ {0,t 1 , ...,t N-1 , T } where the constraints are only satisfied for discrete values over the time-grid. Therefore, the optimization algorithm deals with a set of discrete values corresponding to the continuous constraint values over a time-grid. There are several methods to adapt the grid T in order to get better results [START_REF] Hettich | Semi-infinite programming: theory, methods, and applications[END_REF].

F. The constraints are violated

Using the classical discretization method makes the optimization produces a solution which satisfies the constraints only for the discrete instants {t 0 ,t 1 ...,t N-1 ,t N }, but does not ensure them elsewhere. In order to highlight how such an approach may be hazardous for the robot integrity and balance, we used this method to plan a one-step motion for the HOAP-3 humanoid robot. We considered the constraint (12) and used the motion duration as the objective function and solve this optimization problem thanks to the algorithm C-FSQP [START_REF] Lawrence | User's Guide for CFSQP Version 2.5: A C Code for Solving (Large Scale) Constrained Nonlinear (Minimax) Optimization Problems, Generating Iterates Satisfying All Inequality Constraints[END_REF]. Figure 2 shows the ZMP time-history for an optimal motion obtained with a 10-time-point grid discretization. Here, in order to keep the balance of the robot, the ZMP must stay within the interval [-0, 04; 0, 068] which corresponds to the robot foot size. This constraint is indeed satisfied for the discrete time instants, actually given to the optimization algorithm. However the continuous function is violated, for instance between t = 0.05s and t = 0.1s. To exhibit how this violation may impact the planned motion, we experimented this motion on the HOAP-3 humanoid robot and found that the robot lost its balance (Cf. Fig. 3). One intuitive way to solve this issue could be to increase the number k of discrete time instants. But this will not ensure the avoidance of any constraint violation [START_REF] Lengagne | Guaranteed computation of constraints for safe path planning[END_REF]. In the next section we show how our safe method can solve this issue.

IV. SAFE MOTION PLANNING

A. Definition

Safe motion planning uses the same optimization softwares as usual motion planning. However, the constraints are computed using the guaranteed discretization introduced in our previous work [START_REF] Lengagne | Guaranteed computation of constraints for safe path planning[END_REF]. The main idea at the core of the guaranteed discretization method is to return the extremum values of the inequality constraint functions g i (X,t) during a time interval g i (X, τ) ≤ 0 ( 14)

With IT = {[t] 1 , [t] 2 ], ..., [t] k-1 , [t] k } where time intervals [t] n = [t n-1 ,t n ]
. Interval analysis will allow to do the practical implementation of this method.

B. Interval analysis

Interval analysis was initially developed to account for the quantification errors introduced by the floating point representation of real numbers with computers and was extended to validated numerics [START_REF] Sunaga | Theory of interval algebra and its application to numerical analysis[END_REF], [START_REF] Moore | Methods and Applications of Interval Analysis[END_REF], [START_REF] Neumaier | Interval methods for systems of equations[END_REF]. A real interval

[a] • [b] = [in f u∈[a],v∈[b] u • v, sup u∈[a],v∈[b] u • v] ( 15
)
Consider m : R n -→ R m ; the range of this function over an interval vector [a] is given by:

m([a]) = {m(u) | u ∈ [a]} (16)
The interval function [m] :

IR n -→ IR m is an inclusion function for m if ∀[a] ∈ IR n , m([a]) ⊆ [m]([a]) (17)
An inclusion function of m can be obtained by replacing each occurrence of a real variable by the corresponding interval and each standard function by its interval counterpart. The resulting function is called the natural inclusion function. The performances of the inclusion function depend on the formal expression of m.

Finally, interval analysis offers a practical but guaranteed mean to compute minimum and maximum values for a function m(t) when t is defined over a given interval [

t]. An upper bound for the maximum value max t∈[t] (m(t)) is given by Sup[m]([t]) and a lower bound for the minimum value min t∈[t] (m(t)) is given by In f [m]([t]).

In practice, the bounds thus derived may be too large because of the wrapping and dependence effects. Still, there are several techniques that can be used to obtain tighter ones by using for instance Taylor series expansion or some global optimization techniques [START_REF] Hansen | Global optimization using interval analysis[END_REF].

Therefore the upper bounds for ( 14) max g i are obtained in an easy and practical way by computing the upper bound of the inclusion function

[g i ] for a time interval [t]: ∀i, ∀[t] ∈ IT Sup[g] i (X, [t]) ≤ 0 (18)

C. Practical application

In this study, we have used a simple way to solve the wrapping and dependence effects issue. In order to compute tighter enclosures for the inclusion functions, and hence a better evaluation of the extremum values of a given function over the time-interval, we have chosen to split the timeintervals.

Splitting is the process that decomposes an interval into N b subintervals, then the extremum of the interval may be taken as the extremum of the N b both sub-intervals. Figure 4 shows how the splitting makes the enclosure tighter. However, splitting may increase computation time. These splitting steps are carried out until the estimated maximum and minimum values are obtained with an acceptable accuracy.

As an illustrative example, let us consider the inclusion function [ZMP]([t]) (Equation 5). In the latter equation, the t variable appears several times. Therefore the formal expression of the function ZMP(t) will suffer from the dependence effect. Thanks to the splitting, we can control the accuracy in the computation of the extremum values.

Let us quantify the accuracy of the enclosures as the relative measure of width of the computed enclosures with regard to the actual one, as follows eps(%) = width(computed) -width(actual) width(actual) × 100 [START_REF] Sunaga | Theory of interval algebra and its application to numerical analysis[END_REF] Here, the actual enclosures are the ones obtained with N b = 2 15 subintervals. Table I shows that the bigger the of sub-intervals (N b ), the better the accuracy and the larger the computation time T c . Therefore the user has to make a compromise between accuracy and computation time. In our case we choose to consider 128 sub-intervals.

D. The constraints are satisfied

We experimented a 7cm step motion generated with our safe motion planning method considering 10 intervals, for the humanoid robot HOAP-3. Figure 6 shows this motion, where the robot keeps its balance contrary to figure 3 where it falls. Figure 5 shows both the ZMP time-history ZMP(t) and the guaranteed interval discretization [ZMP]([t] k ) always remain within the base of support. E. The gradient of constraint functions Some optimization algorithms need the computation of the gradient of the constraint functions with respect to the parameter vector ∂ g(X,t) ∂ X . With a grid discretization (Cf Section III-E) this gradient is computed at the grid instant. Now, the method used in subsection IV-C for solving equation ( 18) subdivides a given time interval

∀t ∈ [0, τ] -0, 04 ≤ ZMP(t) ≤ 0, 068 (20)
[t] k into N b subintervals [t] k = ∪ l=1,...,N b [t] k,l
Then there exists a subinterval [t] k,l max which contains the maximum of g i (X,t) for t ∈ [t] k . Since we do not keep track on when the maximum occurs within [t] k,l max , we choose to approximate the gradient of the maximum of g i (X,t) by the value of the gradient obtained at the middle of the subinterval [t] k,l max . This is summarized in the following equation:

∂ ∂ X Sup[g] i (X, [t] k ) ≈ ∂ ∂ X g i (X, Mid[t] k,l max ) (21)
We can infer from this that the number of sub-intervals N b will influence the accuracy of the computed gradient and hence the efficiency of the (C-FSQP or IPOPT) optimization algorithm used. Note also that the partial derivatives can be computed either via formal methods or automatic differentiation [START_REF] Bendtsen | FADBAD, a flexible c++ package for automatic differentation[END_REF]. In our work, we used the latter.

V. RESULTS

A. Procedure

In this section we compare the usual motion planning method with our safe motion planning method. We want to plan a walking motion with steps of 5 cm. Then we consider a cyclic motion where the initial posture is symmetrical to the final one. We consider the 6 dof model of the HOAP-3 humanoid robot (Subsection. II-A). The joint trajectories are computed with 4 B-splines functions. The vector X contains (4 -1) × 6 + 1 = 19 parameters (Here -1 is due to the cyclic nature of the retained motion, i.e. final position is symmetric to the initial one) . We opt for the algorithm IPOPT [START_REF]Introduction to IPOPT : a tutorial for downloading, installing and using IPOPT[END_REF] and perform two optimizations, the first one for minimizing the motion duration and the second one for minimizing the energy consumption. On one hand we asses the classical motion planning method considering a grid of 15 points. On the other hand we use our safe motion planning method considering 5 intervals splitted into 128 sub-intervals.

We repeated 500 times the numerical procedure with random initial conditions to evaluate the impact of our method regarding the repeatability and the time computation of the algorithm.

B. Repeatability

We define the repeatability of a cost function value the percentage of optimization procedure which produce a solution with this value. A high percentage of repeatability means that the optimization algorithm has good convergence properties. Tables II andIII show that there are some differences of the cost function values between the classical and the safe motion planning method (1.055 = 1.075). These differences are mainly due to the fact that the classical motion planning method produces a solution which may violate some constraints, whereas a slower motion will not.

Moreover, one can see that the percentage of a failure of the optimization algorithm is higher with the classical motion planning method (4%) than with the safe motion planning method (0.4%). In the same time the percentage to have one of the two best cost function values is higher considering a safe motion planning (98.6 %) than the classical one (76.8 %). We can see the same phenomena on Tables IV andV, even if the difference is smaller (97.6% against 95.4%). This is due to the fact that an energy-minimum motion is slower than a time-minimum motion and consequently violates less the dynamic constraint functions. Therefore it appears that a guaranteed discretization ensures a better convergence for the optimization algorithm.

C. Computation time

A study about the computation time of these two methods for a simple 2-dof system was done in [START_REF] Lengagne | Guaranteed computation of constraints for safe path planning[END_REF]. Table VI shows the difference of the computation time between the classical and the safe motion planning method considering two different cost functions: the motion duration and the energy consumption. The computation time of the safe motion planning method is slower than the classical one, and the choice of the cost function has an effect on the computation time. Indeed a time-minimum motion has fast dynamics that cause constraint violation, thus the optimization algorithm will call more often the constraint functions than for a slower motion such as energy minimum motion. That is why safe motion planning method is 10 times slower than classical one for a timeminimum motion and only 4 times slower for energy-minimum motion.

D. Experiments

Figure 7 shows an energy-minimum motion computed with the safe motion planning method. One can see that the safe motion planning method generates a motion which avoids any constraint violation since the robot keeps its balance. To the contrary our new safe motion planning method uses Interval Analysis to compute the maximum of the constraint functions over time-intervals, which avoids any constraint violation. We showed that our method increases the repeatability of the algorithm but the computation time is larger.

Here we addressed the one-dimension time discretization issue but the same approach can be used for N-dimensions space discretization. We illustrated our safe motion planning method with a 2-D model of humanoid robots. We plan to use it with a 3-D model.

Fig. 1 .

 1 Fig. 1. 2-D model of the humanoid robot for a one-step motion in the sagital plane

Fig. 2 .

 2 Fig. 2. Representation of the ZMP, for a motion optimized with a 10-timepoints discretization

Fig. 3 .

 3 Fig.3. A 10-time-points discretization time-minimum motion optimization using C-FSQP[START_REF] Lawrence | User's Guide for CFSQP Version 2.5: A C Code for Solving (Large Scale) Constrained Nonlinear (Minimax) Optimization Problems, Generating Iterates Satisfying All Inequality Constraints[END_REF]

  [t] = [t,t] to the optimization algorithm instead of the values taken over the time grid. Then, equation (10) is replaced by: ∀i, ∀[t] ∈ IT max ∀τ∈[t]

  [a] = [a; ā] is a connected and closed subset of R. With a = In f ([a]) and ā = Sup([a]). The set of all real intervals of R is denoted by IR. Real arithmetic operations are extended to intervals. Consider an operator • ∈ {+, -, * , ÷} and [a] and [b] two intervals. Then:

Fig. 4 .

 4 Fig. 4. The bounded ZMP over 6 time-intervals

Fig. 5 .

 5 Fig. 5. Representation of the ZMP, for a motion optimized with 10 intervals

Fig. 7 .

 7 Fig. 7. Energy-minimum 5cm step motion

TABLE V REPEATABILITY

 V FOR AN ENERGY-MINIMUM SAFE MOTION PLANNING