Frédérique Charles
email: charles@ann.jussieu.fr

Bruno Després
email: despres@ann.jussieu.fr

Benoît Perthame
email: benoit.perthame@upmc.fr

Rémi Sentis
email: remi.sentis@cea.fr

Nonlinear stability of a Vlasov equation for magnetic plasmas

Keywords: Vlasov equations, Plasma physics, Kinetic averaging lemma, Maxwell's equations. Mathematics Subject Classification: 35B35, 35L60, 82D10

The mathematical description of laboratory fusion plasmas produced in Tokamaks is still challenging. Complete models for electrons and ions, as Vlasov-Maxwell systems, are computationally too expensive because they take into account all details and scales of magneto-hydrodynamics. In particular, for most of the relevant studies, the mass electron is negligible and the velocity of material waves is much smaller than the speed of light. Therefore it is useful to understand simplified models. Here we propose and study one of those which keeps both the complexity of the Vlasov equation for ions and the Hall effect in Maxwell's equation. Based on energy dissipation, a fundamental physical property, we show that the model is nonlinear stable and consequently prove existence.

Introduction

To describe the behavior of ions population in hot plasmas, it is very classical to address, at least at theoretical level, a Vlasov equation coupled to a non-linear Poisson equation which defines the electrostatic field. When the magnetic field B is an external datum, for the ion distribution function f (t, x, v) (at position x and velocity v), one addresses the following system where n e = n e (t, x) is the electron density, T e is the mean electron temperature assumed to be constant and λ is the Debye length (a characteristic constant of the plasma). This system is a classical one, indeed the electrostatic field may be approximated by -T e ∇ ln n e (at least when the electron temperature is constant) and the Poisson equation is nothing but the Gauss relation applied to this field (cf. [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF] for example)

The mathematical understanding of this kind of kinetic system has made important progresses with the proof of existence of global weak solutions in the large of the Vlasov-Maxwell system in [START_REF] Diperna | Global weak solutions of Vlasov-Maxwell systems[END_REF] and of the Vlasov-Poisson system [START_REF] Lions | Propagation of moments and regularity for the 3-dimensional Vlasov-Poisson System[END_REF][START_REF] Pfaffelmoser | Global classical solutions of the Vlasov-Poisson system[END_REF]. However, the mathematical description of laboratory fusion plasmas (such as those produced in Tokamaks) with this kind of ion kinetic models is still a major challenge, in particular to deal with a time scale compatible with the evolution equation for the magnetic field. As a matter of fact, complete systems of Vlasov-Maxwell type which account for all scales of electro-dynamics are not relevant since the velocity of ion waves is much smaller than the speed of light: it is well known that at the time scale of the ion population it is convenient to neglect the current of displacement in the Maxwell equations (as in magneto-hydrodynamics models cf. [START_REF] Lütjens | The XTOR code for nonlinear 3D simulations of MHD instabilities in tokamak plasmas[END_REF][START_REF] Lütjens | XTOR-2F: A fully implicit Newton-Krylov solver applied to nonlinear 3D extended MHD in tokamaks[END_REF]), that is to say the electric current J is assumed to satisfy

J = ∇ ∧ B.
The aim of this work is precisely to propose, justify and study a model which couples a kinetic equation for the ions and an evolution equation for the magnetic field (as those used in magneto-hydrodynamics). The unknowns are the ion particle density f (t, x, v), the magnetic field B(t, x) and the electron density n e (t, x) and they satisfy

                   -λ 2 ∆ ln n e = n I -n e , (a)
∂B ∂t -∇ ∧ 1 n e n I u I ∧ B + ∇ ∧ 1 n e J ∧ B + ∇ ∧ (η∇ ∧ B) = 0, (b)
∂f ∂t + v • ∇f + ∂ ∂v • - T e n e ∇n e + J -n I u I n e ∧ B + v ∧ B f = 0, (c)
∇ • B = 0, (d) (1)
where the following notations are used

n I (t, x) = R 3 f (t, x, v)dv (the number density in ions), (2)
n I (t, x)u I (t, x) = R 3 f (t, x, v)vdv (the macroscopic velocity of ions).

(

) 3
and η denotes a strictly positive bounded function corresponding to the plasma resistivity. Notice firstly that the electric field E is given by the relation

n e E = -T e ∇n e -n I u I ∧ B + J ∧ B + n e η∇ ∧ B (4)
which is one of the classical forms of the generalized Ohm law [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF][START_REF] Braginskii | Transport Processes in a Plasma[END_REF][START_REF] Blum | Numerical simulation and optimal control in plasma physics with application to Tokamaks[END_REF][START_REF] Freidberg | P lasma physics and fusion energy[END_REF] (it includes the term J ∧ B related to the Hall effect); so, (1)-b) is exactly the Faraday equation ∂B ∂t + ∇ ∧ E = 0. We have chosen here a rescaling, such that the ion mass is set to 1 and the electron and ion charge are also set to 1. Moreover the scaling of the ion and electron density is such that the characteristic value of n I is equal to 1. So with this scaling, the Debye length λ is defined by λ 2 = T e ε 0 where ε 0 is vacuum dielectric constant and equation (1)-a) reads also as -ε 0 T e ∆ ln n e = n I -n e which is Gauss relation applied to -T e ∇ ln n e , the dominant term in the electric field. In the rest of the paper, we use the following notations

E 0 (t, x) = -T e ∇n e n e + 1 n e (J -n I u I ) ∧ B and F(t, x, v) = E 0 + v ∧ B. ( 5
)
Here F is the kinetic force field which appears in equation (1)-c).

Remark 1. It is fundamental to notice that the Debye length is a small quantity in many situations; it may be proved that in the limit λ → 0 the solution n e to the non-linear Poisson equation satisfies n e → n I and many studies have analyzed this singular limit in different contexts [START_REF] Brezis | Analyse asymptotique de l'équation de Poisson couplée à la relation de Boltzmann. Quasi-neutralité dans les plasmas[END_REF][START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF][START_REF] Han Kwan | Quasineutral limit of the Vlasov-Poisson system with massless electrons[END_REF]. Then, if we make the approximation n e ≃ n I , the Ohm law (4) leads to the classical formula

E + u I ∧ B = -T e ∇ ln n I + 1 n I (∇ ∧ B) ∧ B + η∇ ∧ B
and from (1), we recover formerly the following system

∂B ∂t -∇ ∧ (u I ∧ B) + ∇ ∧ 1 n I (∇ ∧ B) ∧ B) + ∇ ∧ (η∇ ∧ B) = 0, ∂f ∂t + v.∇f + ∂ ∂v (E lim + v ∧ B)f = 0, n I (E lim + u I ∧ B) = -T e ∇n I + (∇ ∧ B) ∧ B.
Despite its apparently simpler form than the system (1), it is not used in physical literature, up to our knowledge. As a matter of fact even if the magnetic field is neglected, the equation

∂f ∂t + v.∇f + ∂ ∂v (E lim + v ∧ B)f = 0, n I E lim = -T e ∇n I
is probably ill-posed from a mathematical point of view (nevertheless [START_REF] Han Kwan | Quasineutral limit of the Vlasov-Poisson system with massless electrons[END_REF] has noticed that using some arguments of the paper [START_REF] Mouhot | On Landau damping[END_REF], one may find a weak solution on a small time interval). Moreover it would be hard to give a meaning of the quantity u I = R 3 f (v)vdv/ R 3 f (v)dv which appears in the magnetic equation (it is hardly possible to derive a bound of n I away from below). Notice that this kind of model has been addressed in [START_REF] Ghendrih | D erivation of a gyrokinetic model. Existence and uniqueness of specific stationary solution[END_REF] and [START_REF] Acheritogaray | Kinetic formulation and global existence for the Hall-Magneto-hydrodynamics system[END_REF] but the framework is a little different since a velocity diffusion term (of Fokker-Planck type) is accountered for in the kinetic equation.

Back to (1), from a mathematical point of view, one can show that n e is naturally bounded away from zero due to the nonlinear elliptic equation (1)-a and the quantity nI ne u I = n -1 e R 3 f (v)vdv may be defined in some Lebesgue space.

Our study is aiming at the analysis of model ( 1) in a bounded domain; therefore we have to specify the boundary and initial conditions.

Boundary conditions

Let Ω be a bounded domain of R 3 of class C 1,1 . We consider the model (1) for v ∈ R 3 , x ∈ Ω and t ≥ 0, with the following boundary conditions

     n x • ∇n e (t, x) = 0 x ∈ ∂Ω, (a)
n x ∧ B(t, x) = 0 x ∈ ∂Ω, (b) f (x, v -2(v • n x )n x ) = f (x, v), v ∈ R 3 , x ∈ ∂Ω. (c) (6)
The zero flux boundary conditions (6)-a) on electrons enforces in (1)-a) the global neutrality of the plasma

Ω n e dx = Ω n i dx.
The condition (b) on the magnetic field is the simplest one and a more realistic condition is discussed in section 4. Our result will be still true with more realistic boundary condition, of impedance type for example, or also a non homogeneous boundary condition such as

n x ∧ B(t, x) = n x ∧ B imp , x ∈ ∂Ω,
which is much more relevant in the context of confined plasmas in Tokamaks. Condition (6)-c) is the so-called specular reflection; it implies the no-slip boundary condition

u I (t, x) • n x = 0, x ∈ ∂Ω. (7)
Our main result is the following.

Theorem 1 (Nonlinear stability). Assume [START_REF] Freidberg | P lasma physics and fusion energy[END_REF] and that the initial datum satisfy (11)- [START_REF] Evans | Partial Differential Equations[END_REF]. From any bounded family of solutions f ε , B ε , n ε e to (1), [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF] we can extract a subsequence that converges to a weak solution with finite energy. Moreover, we have for all T > 0,

B ε → B ∈ L ∞ 0, T ; L 2 (Ω) ∩ L 2 0, T ; L q (Ω)
strongly, q < 6,

f ε → f ∈ L ∞ 0, T ; L 1 ∩ L ∞ (Ω × R 3 )) weakly, J ε = ∇ ∧ B ε → J in L 2 ((0, T ) × Ω) weakly,
n ε e → n e in L q ((0, T ) × Ω), strongly and for some constants 3 strongly, for all q such that 1 ≤ q < 5/4.

0 < K -≤ n e ≤ K + , n ε I u ε I → n I u I ∈ (L q (0, T × Ω))
The first tool towards a mathematical analysis is to state the energy balance for the full system; this is performed in section 2. It indicates that the magnetic resistivity in equation (1)-a) is useful (if not necessary) so as to control the current J and subsequently the term J ∧ B describing to the Hall effect. Also the classical tools for elliptic equation with Neumann boundary conditions will lead to control n e away from below. But more technical ingredients (namely a kinetic averaging lemma in a L q space) are used in the proof in order to pass to the limit in the nonlinearities; this is the purpose of section 3 and of several appendices. We come back on more realistic boundary conditions in section 4 so as to include confinement. We also detail a constructive splitting procedure which can be used to design an approximate solution with the same energy law.

Preliminaries

The first and major ingredient in our approach to nonlinear stability is the energy balance. We state it here together with considerations on the physics sustaining the model.

The energy balance

We now turn to the study of energy dissipation for the system (1) with homogeneous boundary conditions [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF]. We introduce the ions kinetic energy and the magnetic energy

E I (t) = 1 2 Ω R 3 f (t, x, v)|v| 2 dvdx, E m (t) = 1 2 Ω |B(t,
d dt E I + E m + λ 2 2 Ω |∇ (ln n e )| 2 dx + Ω (n e ln n e -n e + 1) dx = - Ω η ∇ ∧ B 2 dx. (9)
Proof. Notice that, as usual, the free divergence condition (1)-(d) is not needed for the energy identity. We first consider the ion kinetic energy. Using the specular reflection condition, we compute

d dt E I = Ω R 3 F(t, x, v) • vdv = Ω R 3 1 n e (∇ ∧ B) ∧ B + (v - n I n e u I ) ∧ B -T e ∇ (ln n e ) f • vdvdx = Ω n I n e J ∧ B -n I T e ∇ (ln n e ) • u I dx.
Now, we turn to the magnetic energy and first recall the definition of E in (4). Thanks to a classical tensorial identity, we find

d dt E m = Ω ∂B ∂t • Bdx = - Ω ∇ ∧ E • Bdx = - Ω [E • ∇ ∧ B + ∇ • (E ∧ B)]dx = - Ω - n I n e u I ∧ B + η∇ ∧ B • (∇ ∧ B) dx - ∂Ω n x • (E ∧ B)dx.
According to the boundary condition ( 6)-(b), we have

n x • (E ∧ B) = -E • (n x ∧B) = 0, then the boundary integral is zero. Now, thanks to identity J • (u I ∧ B) = -u I • (J ∧ B), we get d dt E m = - Ω n I n e u I • (J ∧ B)dx - R 3 η ∇ ∧ B 2 dx.
At this stage, we have obtained

d dt [E i + E m ] = Ω [(J ∧ B) -n I ∇ (ln n e )] • u I dx - Ω u I • (J ∧ B)dx - Ω η ∇ ∧ B 2 dx = -Ω n I ∇ (ln n e ) • u I dx -Ω η ∇ ∧ B 2 dx. (10)
To continue, we recall the conservation of mass on n

I ∂n I ∂t + ∇ • (n I u I ) = 0.
Then, we may use the boundary condition ( 7) to obtain

- Ω n I ∇ (ln n e ) • u I dx = Ω (ln n e )∇ • (n I u I ) dx = - Ω (ln n e ) ∂n I ∂t dx = - Ω (ln n e ) ∂n e ∂t -λ 2 ∂ ∂t ∆ (ln n e ) dx = - d dt Ω (n e (ln n e ) -n e ) dx -λ 2 d dt Ω |∇(ln n e )| 2 2 dx
where we have used equation ( 1)-(a) once differentiated in time.

Altogether, we obtain the relation (9).

Physical considerations

Let first focus on the momentum balance. To find it, we multiply the kinetic equation ( 1)-c) by v, integrate over R 3 and define as usual [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF] the ion pressure tensor

P I through ∇ • (v ⊗ vf (v))dv = ∇ • (n I u I ⊗ u I ) + ∇ • P I .
Then, we get

∂(n I u I ) ∂t + ∇ • (n I u I ⊗ u I ) + ∇ • P I + T e ∇n e n e n I + n I n e n I u I ∧ B - n I n e J ∧ B -n I u I ∧ B = 0,
that is to say

∂(n I u I ) ∂t + ∇ • (n I u I ⊗ u I ) + ∇ • P I + T e ∇n e -J ∧ B = (n e -n I ) n e [T e ∇n e + (n I u I -J) ∧ B] .
On the left hand side one recognizes the classical momentum equation which appears in the MHD modeling with the magnetic pressure tensor. Let us stress that with the electric field E 0 introduced in (5) and according to (1)-a), the r.h.s. term reads as follows

ε 0 ∇ • E 0 + ∇ • n I u I -J n e ∧ B E 0 = ε 0 ∇• 1 2 |E 0 | 2 -E 0 ⊗ E 0 + λ 2 E 0 ∇ • n I u I -J n e T e ∧ B .
Thus we may see that this equation is in a conservative form only when the Debye length vanishes. This is due to the asymptotic expansion motivating ( 1); we have chosen to keep energy dissipation and smoothness for n e to the expense of momentum balance. In order to ensure exact conservation momentum, a possible route is to use a modified Poisson equation as follows

-λ 2 ∆ ln n e = n I -n e -λ 2 ∇ • J -n I u I n e T e ∧ B
which corresponds to an approximation at the same order but which looses nice properties on n e .

Remark 2 (Electron momentum balance). We see that E 0 is equal to the electric field E up to the resistive part η∇ ∧ B. Let us stress that (∇ ∧ B)-n I u I = J -n I u I is the electron contribution to the electric current and (4) may be interpreted as the electron momentum balance equation when using the massless electron approximation.

Remark 3 (Friction). In order to account for the friction between ions and electrons in the kinetic equation, we could add a friction term to equation (1)-c) and arrive to

∂f ∂t + v.∇f + ∂ ∂v .(Ff ) = - ∂ ∂v (ηJf )
(recall that the coefficient η is related to the collision frequency of electrons against ions). Also, for the sake of compatibility, instead of equation ( 1)-a), the Poisson equation would read

λ 2 -∆ ln n e + 1 T e ∇(ηJ) = n I -n e .
3 Proof of the stability theorem

We give a family of initial conditions that satisfy

f ε (0, x, v) = f in,ε (x, v) ≥ 0, f in,ε is bounded in L 1 ∩ L ∞ (Ω × R 3 ), (11)
B ε (0, x) = B in,ε (x), ∇ x • B in,ε = 0. ( 12
)
From the Lemma 1 in Appendix A and the explanations in Appendix B we see that if E I (0) is uniformly bounded then the free energy at initial time

Ω |∇ x (ln n ε e (0, x))| 2 dx + Ω n ε e (0, x) ln n ε e (0,
x)dx is also uniformly bounded. Then we only need to assume that sup

0<ε≤1 [E I (0) + E m (0)] < ∞ (13)
to have a total energy uniformly bounded for all time. Moreover we assume a control on resistivity as

0 < η min ≤ η ∈ L ∞ (Ω). (14)
Associated with such initial datum, we consider a sequence of strong solutions f ε , B ε , and n ε e of (1) with boundary conditions (6).

A priori bounds

Before we begin the proof of Theorem 1, we recall several general a priori estimates that are used. The energy bound [START_REF] Evans | Partial Differential Equations[END_REF] is of course at the heart of our analysis because it allows us to deduce directly from the estimate (9) the Proposition 2 (Bounds derived from energy). Under the assumptions (11)-( 14), the sequences f ε , B ε , and

n ε e satisfy sup 0<ε≤1 sup t∈[0,∞] R 3 Ω f ε i (t, x, v)|v| 2 dvdx < ∞, (15)
sup 0<ε≤1 sup t∈[0,∞] Ω |B ε (t, x)| 2 dx < ∞, (16)
sup 0<ε≤1 sup t∈[0,∞] Ω |∇ x (ln |n ε e |)| 2 + n ε e ln |n ε e | + n ε e dx < ∞. (17)
Indeed, for [START_REF] Golse | Regularity of the moments of the solution of a transport equation[END_REF], we use the fact that the map x ∈ R → x ln(x) -x is bound from below and thus all the terms in ( 9) are bounded by the initial energy plus a constant. This provides a control on the kinetic energy for the ion distribution that is essential to obtain further a priori bounds Proposition 3 (Estimates for the kinetic densities). For strong solutions one has

f ε ∈ L ∞ 0, ∞; L 1 ∩ L ∞ (Ω × R 3 ) and f ε (t, •, •) p = f in,ε (•, •) p ∀t ≥ 0, ∀p ∈ [1, ∞], n ε I (t, •) 5/3 ≤ C f in,ε (•, •) 2/5 ∞ Ω R 3 f (t, x, v)|v| 2 dvdx 3/5 , n ε I u ε I (t, •) 5/4 ≤ C f in,ε (•, •) 1/5 ∞ Ω R 3 f (t, x, v)|v| 2 dvdx 4/5 . ( 18
)
6

The bounds of f ε follow immediately from the observation that div v F = 0. The bounds on n ε I and n ε I u ε I follow from interpolation inequalities that we recall in the Appendix A.

Next we use the L 5/3 integrability of n I to obtain Proposition 4 (Estimates on the electron density). The electron density satisfies for some constants K + > 0, K -> 0 (depending on the bounds stated as now)

0 < K -≤ n e (t, x) ≤ K + . ( 19
)
This is an easy consequence of elliptic regularity and we give a proof in the Appendix B, it is just a combination of Lemma 1 and Lemma 2.

We now come to the magnetic field. We have the following Proposition :

Proposition 5 (Estimate on the magnetic field). The family of magnetic fields satisfies for a uniform constant C

sup 0<ε≤1 ∞ 0 Ω ∇ x ∧ B ε 2 dxdt ≤ C, sup 0<ε≤1 B ε L 2 0,∞;L 6 (Ω) ≤ C. ( 20
)
Proof. The control in L 2 (0, ∞; H curl (Ω)), follows from the energy dissipation in Proposition 1. Then thanks to Theorem 2, the family B ε is bounded in L 2 (0, ∞; H 1 (Ω) 3 ), and then, thanks to the Sobolev injections, the family B ε is bounded in L 2 (0, ∞; L 6 (Ω) 3 . We refer to Appendix C for precise statements and references.

Our last task is to show that we have enough bounds to define the electric and force fields. The four terms that compose the electric field

E ε = -T e ∇ ln n ε e - n ε I u ε I n ε e ∧ B ε + J ε ∧ B ε n ε e + η∇ ∧ B ε
are respectively uniformly bounded in the spaces

• T e ∇ ln n ε e ∈ L ∞ t (L 2
x ) (energy inequality),

• n ε I u ε I n ε e ∧B ε ∈ L 2 t (L 30/29
x ); this follows from ( 18), ( 19) and ( 20) and the fact that, according to Hölder inequality, we get the bound

T 0 [ AB 30/29 ] 2 dt ≤ T 0 [ A 5/4 B 6 ] 2 dt ≤ sup t A 2 5/4 T 0 [ B 6 ] 2 dt. • J ε ∧ B ε n ε e ∈ L 1 t (L 3/2 x ) ∩ L 2 t (L 1
x ) from the Hölder inequality with [START_REF] Lions | Propagation of moments and regularity for the 3-dimensional Vlasov-Poisson System[END_REF] and the a priori estimates in Proposition 5. By interpolation, we also find that

J ε ∧ B ε n ε e ∈ L r t (L p x ), 1 ≤ r ≤ 2, 1 p = 4 3 - 2 3r .
In particular

J ε ∧ B ε n ε e ∈ L 20/11 t (L 30/29 x ). • η∇ ∧ B ε ∈ L 2 t,x (energy dissipation).
Therefore it is well defined and the four terms will have a weak limit in Lebesgue spaces.

For the kinetic force field,

F ε (t, x, v) := -T e ∇ ln n ε e + v - n ε I n ε e u ε I ∧ B ε + J ε ∧ B ε n ε e ,
the same integrability, locally uniformly in v, holds true because these are the same terms and thus they are bounded as for the electric field at the exception of v ∧ B ε which is well defined thanks to (16).

Space-time compactness

The bounds mentioned before allow us to extract subsequences that converge weakly and our task is now to prove several strong convergence results.

We first observe that, from Proposition 5 and ( 16), we may extract a subsequence, still denoted by B ε such that

B ε -→ ε→0 B strongly in L p (0, T ; L q (Ω) 3 , 1 ≤ p < 2, 1 ≤ q < 6, and 1 ≤ p < ∞, 1 ≤ q < 2. ( 21
)
This is because (see Appendix C) thanks to our a priori estimates and Theorem 2, the field {B ε } is bounded in L 1 (0, T, H 1 (Ω) 3 ), and then, thanks to Rellich Theorem, {B ε (t, •)} is relatively compact in L 6 (Ω) 3 for all t ∈]0, T [).

Compactness in time then follows from the Simon-Lions-Aubin lemma (see [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]) because, as proved at the end of section (3.1), we control L r t (L p x ) integrability of the electric field and

∂B ε ∂t = ∇ ∧ E ε .
According to the kinetic averaging lemma recalled in Appendix D, we know that after extraction of a subsequence, for every test function ψ = ψ(v) compactly supported, there exists a weak limit f * of f ε such that

ψf ε → ψf *
a.e. and ψvf ε → ψvf * a.e.

Now, since

T 0 Ω | ψf ε | 5/3 dxdt and T 0 Ω
| ψvf ε | 5/4 dxdt are uniformly bounded, according to classical interpolation arguments we have for given exponents q and r (with 1 < q < 5/3 and 1 < r < 5/4)

T 0 Ω | ψf ε -ψf * | q dxdt → 0, and T 0 Ω | ψvf ε -ψvf * | r dxdt → 0.
Therefore we get

n ε I -→ ε→0 n I , strongly in L q (0, T ; Ω), (22)
n I u ε I -→ ε→0 n I u I strongly in L r (0, T ; Ω) 3 . ( 23
)
Next for the electronic density we may extract a subsequence, still denoted n ε e such that,

n ε e -→
ε→0 n e a.e. and strongly in L q (0, T ; Ω) 3 .

This is a consequence of the bounds ( 17) and ( 19) and the strong convergence [START_REF] Mouhot | On Landau damping[END_REF].

Passing to the limit

We recall that a weak (distributional) solution to (1) is defined by testing against smooth test functions Ψ(t, x), Φ(t, x, v) respectively for (1)-b) and (1)-c) and satisfying the corresponding boundary conditions as in [START_REF] Cercignani | The mathematical theory of dilute gases[END_REF].

After integration by parts, we obtain respectively for (1)-b) and (1)-c), the definitions

- (0,T )×Ω ∂Ψ(t, x) ∂t B ε + E ε (t, x) • ∇ ∧ Ψ(t, x) dtdx = Ω B in,ε (x, v)Ψ(0, x), (25)
and

- (0,T )×Ω×R 3 ∂Φ(t, x, v) ∂t +v • ∇ x Φ(t, x, v) + F ε (t, x, v) • ∇ v Φ(t, x, v) f ε (t, x, v)dtdxdv = Ω×R 3 f in,ε (x, v)Φ(0, x, v). (26)
The elliptic equation (1)-a) is more standard and does not yield difficulties, hence we do not consider it here. In order to conclude the proof of Theorem 1, our purpose is to show that the limit as ε → 0 of the various unknowns B ε , E ε (t, x), F ε and f ε still satisfy these equalities.

Additionally to the strong convergence results stated in ( 21)-( 24), we are know other weak convergences for functions of interest. Firstly, we have

J ε -⇀ ε→0 J in L 2 ((0, T ) × Ω) 3 .
This is enough to pass to the limit weakly in the electric fields because the nonlinear terms are always formed of either strongly convergent terms or J ε multiplied by a term that converges strongly. As a conclusion, we have

E ε -⇀ ε→0
E in L 20/11 0, T ; L 30/29 (Ω) 3 .

The same applies to the force field. However for later purposes, it is better to use the splitting

F ε = E 0,ε + v ∧ B ε ,
and to notice that, as before, E 0,ε -⇀ ε→0 E 0 in L 20/11 0, T ; L 30/29 (Ω)

3
while the other term converges strongly.

These observations allow us to pass to the weak limit in equations (1)-a) and ( 1)-b), on the electronic density and on the magnetic field. For passing to the weak limit in the third equation (1)-c), let us introduce two smooth test functions χ and ϕ with compact support respectively in [0, T ) × Ω and in R 3 ; then we have to deal with two terms. The first one which reads as

T 0 Ω R 3 B 0,ε ∧ v • ∂ϕ ∂v (v)χ(t, x)f ε (t, x, v)dtdxdv,
can be treated as usual because the magnetic field converges strongly as stated earlier. The second one reads as follows

T 0 Ω R 3 E 0,ε • ∂ϕ ∂v (v)f ε (t, x)χ(t, x)dtdxdv = T 0 Ω χ(t, x)E 0,ε • ψf ε (t, x)dtdx
with ψ i = ∂ϕ ∂vi ; we need to use the kinetic averaging lemma which is stated in appendix D. It states that, after extraction, velocity averages converge strongly

χ ψ i f ε (t, x) → χ ψ i f * in L q ([0, T ] × Ω)
for a given value of the index q, with 1 ≤ q < p. Therefore they converge almost everywhere and this is enough to pass to the weak-strong limit in the term

T 0 Ω χE 0,ε • ψf ε dtdx.
This concludes the proof of the Theorem 1.

Non homogeneous magnetic boundary condition

In view of applications to the physics of confined plasmas, it is worthwhile considering, instead of ( 6)-(b), a non homogeneous boundary condition for the magnetic field

n x ∧ B(t, x) = n x ∧ B imp , x ∈ ∂Ω, (27)
The imposed magnetic field B imp (x) is given and can be a function of time as well eventhough we restrict our analysis to space dependency only. In order to introduce this boundary condition in the problem, we assume that B imp can be defined globally as a smooth function in the closure of Ω. It turns out that a convenient regularity assumption is

B imp 1,∞ ≤ C. (28)
We will write B = B imp + B pert where B pert is the perturbation. Using these notations the system (1)-(a,b,c) is rewritten as

                       -λ 2 ∆ ln n e = n I -n e , (a)
∂B pert ∂t -∇ ∧ 1 n e n I u I ∧ (B imp + B pert ) + ∇ ∧ 1 n e (J imp + J pert ) ∧ (B imp + B pert ) +∇ ∧ (η (J imp + J pert )) = 0, (b) ∂f ∂t + v.∇f + ∂ ∂v (- T e n e ∇n e + (J imp + J pert ) -n I u I n e ∧ (B imp + B pert )) + v ∧ (B imp + B pert ) f = 0, (c).
(29) where the total current is J = J imp + J pert with

J imp = ∇ ∧ B imp and J pert = ∇ ∧ B pert .
We now use the boundary conditions

     n x • ∇n e (t, x) = 0 x ∈ ∂Ω, (a)
n x ∧ B(t, x) pert = 0 x ∈ ∂Ω, (b) f (x, v -2(v • n x )n x ) = f (x, v), v ∈ R 3 , x ∈ ∂Ω. (c) (30)
Let us define the perturbed magnetic energy and the total perturbed energy

E pert m (t) = 1 2 Ω |B pert (t, x)| 2 dx, E pert tot = E I + E pert m + λ 2 2 Ω |∇ x (ln n e )| 2 dx + Ω (n e ln n e -n e + 1) dx ≥ 0.
The energy balance is modified by the imposed magnetic field and we have Proposition 6. Classical solutions to (29)-(30) satisfy the perturbed energy dissipation relation

d dt E pert tot = - Ω η J pert 2 dx + S
where the source is

S = - Ω η J imp • J pert dx - Ω 1 n e J imp ∧ (B imp + B pert ) • J pert dx + Ω 1 n e J imp ∧ (B imp + B pert ) • n I u I dx.
Proof. Performing the same manipulations as in the proof of the energy identity [START_REF] Dautray | Mathematical analysis and numerical methods for sciences and technology[END_REF], we obtain

d dt E I = Ω n I n e (J imp + J pert ) ∧ (B imp + B pert ) -n I T e ∇ (ln n e ) • u I dx.
Taking the scalar product of the magnetic equation against B pert and using the homogeneous boundary condition (30)-(b), we obtain

d dt E pert m = - Ω n I n e u I • (J pert ∧ (B imp + B pert )) dx - Ω 1 n e (J imp + J pert ) ∧ (B imp + B pert ) • J pert dx - R 3 η (J imp + J pert ) • J pert dx = - Ω n I n e u I • (J pert ∧ (B imp + B pert )) dx - Ω 1 n e J imp + ∧ (B imp + B pert ) • J pert dx - R 3 η (J imp + J pert ) • J pert dx
Therefore one gets

d dt E I + E pert tot = - Ω T e ∇ (ln n e ) • (n I u I )dx - Ω η J pert 2 dx + S
which is very similar to [START_REF] Després | Reduced resistive MHD in Tokamaks with general density[END_REF]. The rest of the proof is unchanged.

Proposition 7. There exists a constant K > 0 depending only on the initial data and on the constant in [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF] such that the perturbed energy is bounded for all time t < T

⋆ = log 1 + C Jimp ∞ .
Remark 4. Our analysis of weak stability can easily be extended to this non-homogeneous boundary condition for t < T ⋆ . This estimate expresses the interest of a good control on the imposed current. Indeed the smaller is J imp ∞ , the greater T ⋆ is.

Proof. (of proposition 7). We observe that the new terms in S are proportional to J imp . Therefore two cases occur.

First case: J imp = 0. This idealized case might be encountered in Tokamaks: for example if B imp = F ∇θ where θ is the toroidal angle and F is a constant. Physically it corresponds to an imposed exterior magnetic with vanishing current [START_REF] Després | Reduced resistive MHD in Tokamaks with general density[END_REF]. The magnetic lines of B imp form a ring. In this case S = 0, so the perturbed energy dissipation relation has the same form as [START_REF] Dautray | Mathematical analysis and numerical methods for sciences and technology[END_REF]. It turns out that the perturbed energy is bounded for all times, which indeed corresponds to the claim since J imp = 0.

Second case: J imp = 0. This situation is more realistic in Tokamaks, since the magnetic lines form an helix or at least must be close to helicoidal geometry [START_REF] Després | Reduced resistive MHD in Tokamaks with general density[END_REF]. This is the general case.

Let σ > 0 be a given positive number. The energy identity writes also

d dt e -σt E pert tot = Re -σt , R = -σE pert tot - Ω η J pert 2 dx + S. ( 31
)
Our goal is to control the source term S in R as much as possible by -σE pert tot -η J pert 2 and to control the remaining part with a Gronwall technique.

Using the assumption ( 28), the fact that η is constant and the Hölder inequality between the conjugated spaces L 5 (Ω) and L 5 4 (Ω), we can write

|S| ≤ α 1 J pert 2 + α 2 J pert 2 B pert 2 + α 2 B imp 5 n I u I 5 4 + α 2 B pert 5 n I u I 5 4 J imp ∞ with α 1 = η and α 2 = n -1 e ∞ .
Let ǫ > 0 be an arbitrary real number. Then

α 1 J pert 2 + α 2 J pert 2 B pert 2 ≤ α 1 2ǫ + ǫ α 1 + α 2 2 J pert 2 2 + α 2 2ǫ B pert 2 2 .
The identity (34), written as n

I u I 5 4 ≤ α 4 E 4 5
I , implies successively the controls

B imp 5 n I u I 5 4 ≤ α 3 E 4 5 I B pert 5 n I u I 5 4 ≤ α 5 ( B pert 2 + J pert 2 ) E 4 5 I ≤ α 6 E pert m 1 2 E 4 5 I + ǫ α 5 2 J pert 2 2 + α 5 2ǫ E 8 10
I .

So the right hand side in (31) is bounded by

R ≤ β 1 + β 2 E 4 5 I + β 3 E pert m 1 2 E 4 5

I

for some constants β 1,2,3 . Since the electronic density is bounded and

0 ≤ E I + E pert m ≤ E pert tot by construction, we also have R ≤ γ 1 + γ 2 E pert tot 4 5 + γ 3 E pert tot 13 10
for some constants γ 1,2,3 which do not depend on time. Since y for some constants δ 1,2 which do not depend on time.

We see that the right hand side is more than linear with respect to E pert tot due to the power 13 10 . As a consequence this inequality cannot prove that e -σt E pert tot or E pert tot is bounded for all time. Next we wish to obtain a evaluation of the time of existence with respect to J imp ∞ . We set u(t) = e -σt E pert tot . One can simplify the inequality as

u ′ (t) ≤ δ 1 + δ 2 u 13 10
J imp ∞ e 3σ 10 t .

Rescaling of the time variable as dτ = J imp ∞ e 3σ 10 t dt, that is

τ = 10 3 J imp ∞ e 3σ 10 t -1 , yields the inequality d dτ u ≤ δ 1 + δ 2 u 13 10 . It is finally convenient to define v = δ 1 + δ 2 u 13 10 so that v ′ (t) = δ 2 13 10 u 3 10 u ′ (t) which yields v ′ (t) ≤ δ 3 v 3 13 v = δ 3 v 16 13 . Therefore -d dτ v -3 13 ≤ δ 4 = 3 13 δ 3 . which implies v -3 13 (0) -v -3 13 (τ ) ≤ δ 4 τ . It yields v 3 13 (τ ) ≤ 1 v 3 13 (0) -δ 4 τ which is valid for τ < τ ⋆ = v 3 13 (0)
δ4 . Going back to the time variable t, the solution is defined for t < T ⋆ where 10 3

J imp ∞ e 3σ 10 T ⋆ -1 = τ ⋆ .
The proof is complete.

Construction of an approximate solution

In order to complete our theory, we now detail how to use the so-called splitting strategy, which is a constructive method, for the design of an approximate solution to the system

               -λ 2 ∆ln n e = n I -n e , (a)
∂B ∂t -∇∧ 1 ne n I u I ∧ B + ∇ ∧ ( 1 ne J ∧ B)+∇ ∧ (η∇ ∧ B) = 0, (b)
∂f ∂t + v.∇f + ∂ ∂v (- T e n e ∇n e + J -n I u I n e ∧ B) + v ∧ B f = 0, (c). (32)
The idea is clearly inspired from numerical methods. It consists of a convenient splitting strategy à la Strang, together with the linearization and freezing of certain coefficients à la Temam. The main point is to decompose the total system in simpler parts which are conceptually easier to solve or easier to analyse, preserving at the same time the decay of the energy identity

E tot = 1 2 Ω R 3 f (t, x, v)|v| 2 dvdx + 1 2 Ω |B(t, x)| 2 dx + λ 2 2 Ω |∇ x (ln n e )| 2 dx + Ω
(n e ln n e -n e + 1) dx.

Let ∆t > 0 be a time step which is ultimately destinated to tend to zero. We consider that

f (t k ) and B(t k )
are known at the beginning of the time step t k = k∆t. We restrict the presentation to the core of the method. This constructive method also provides additional insights into the mathematical structure of the model.

Vlasov-Poisson

One first solves during the time step ∆t

               -λ 2 ∆ln n e = n I -n e , (a)
∂B ∂t = 0, (b) ∂f ∂t + v.∇f + ∂ ∂v - T e n e ∇n e f = 0. (c)
This is a non linear Vlasov-Poisson equation which can be considered as standard even if we know very little mathematical literature about it. It is easy to show that regular solutions preserve the energy. This procedure defines a new solution f ⋆ (t k + ∆t) and B ⋆ (t k + ∆t) = B(t k ).

Magnetic part, first stage

For convenience we split the magnetic part of the equations in two stages, the first one which is fundamental, and the second which deals with less involved terms. The first stage writes

         ∂B ∂t -∇∧ 1 ne n I u I ∧ B frozen + ∇ ∧ ( 1 ne J ∧ B frozen )+∇ ∧ (η∇ ∧ B) = 0, (b) ∂f ∂t + ∂ ∂v J n e ∧ B frozen ) f = 0, (c).
The initial data is provided by the previous step of the algorithm

f (t k ) = f ⋆ (t k + ∆t) et B (t k ) = B ⋆ (t k + ∆t).
The frozen magnetic field is

B frozen = B ⋆ (t k + ∆t).
This frozen field is constant in time during the whole time step. This trick was first introduced by Temam in the context of trilinear forms and Navier-Stokes equations for magnetic equations [START_REF] Temam | Remarks on a free boundary value problem arising in plasma physics[END_REF][START_REF] Temam | N avier-Stokes Equations, Theory and Numerical Analysis[END_REF], see also [START_REF] Gerbeau | Mathematical methods for the magnetohydrodynamics of liquid metals[END_REF]. One notices that equation b) is now a linear one, even if equation c) is still formerly non linear because it has a n e dependence. However an explicit procedure allows to compute the solution. Indeed solutions of equation c) are such that

∂ t n I = 0.
It means that n I and n e are frozen quantities

n I = n I frozen et n e = n e frozen .
One has Once the magnetic field is computed, we can report the current J in equation c) which is now easily solved with the method of characteristics. It is immediate that, due to the resistive operator, the energy decreases during this step. Since n i is constant, the electronic energy is constant.

∂ t n I u I = J ∧ d, d = n I frozen B frozen
The solution at the end of this stage is referred to as

f (t k + ∆t) and B (t k + ∆t) = B(t k ).

Magnetic part, second stage

It remains to solve

         ∂B ∂t = 0, (b) ∂f ∂t + ∂ ∂v -n I u I n e ∧ B) + v ∧ B f = 0, (c), with prescribed initial data f • (t k ) = f (t k + ∆t) and B • (t k ) = B (t k + ∆t).
Since the magnetic field is frozen B = B frozen , equation c) is greatly simplified. Once again the ionic density n I and the electronic density n e are frozen n I = n I frozen and n e = n e frozen .

A consequence is

∂ t n I u I = n I u I ∧ d, d = - n I frozen n e frozen + 1 B frozen .
The solution of this linear equation is immediate. Therefore n I u I is known. And finally the method of characteristics can be used to solve c). The energy is preserved during this second magnetic stage.

Iterations

The previous procedure allows us to design an approximate solution

f (t k + ∆t) = f • (t k + ∆t) and B(t k + ∆t) = B • (t k + ∆t)
one time step after the other. With this procedure the total energy decreases and we can apply our stability analysis for proving existence.

A Control on moments of f I Several type of controls on velocity moments of f I are available, see [START_REF] Perthame | Mathematical Tools for Kinetic Equations[END_REF]. Here we recall one of the most fundamental control in L p spaces based on the kinetic energy.

Lemma 1. Let f ∈ L ∞ t,x,v ((0, T ) × Ω × R 3 ) ∩ L ∞ t (0, T ; L 1 x,v (Ω × R 3 , |v| 2 dxdv). Define n I by (2) and n I u i by (3). Then n I ∈ L ∞ t 0, T ; L 5/3 x (Ω) , n I u I ∈ L ∞ t 0, T ; L 5/4
x (Ω) and we have for all t ∈ [0, T ],

n I (t, •) 5/3 ≤ C f (t, •, •) ∞ 2/5 Ω R 3 f (t, x, v)|v| 2 dvdx 3/5 , (33) n
I (t, •)u I 5/4 ≤ C ′ f (t, •, •) ∞ 1/5 Ω R 3 f (t, x, v)|v| 2 dvdx 4/5 . (34)
Proof. We only recall the proof of the result on n I . Let R > 0. We have

n I (t, x) = |v|≤R f (t, x, v)dv + |v|≥R f (t, x, v)dv ≤ CR 3 f (t, •, •) ∞ + 1 R 2 |v|≤R f (t, x, v)|v| 2 dv.
Then by minimization over R we get

n I (t, x) ≤ C f (t, •, •) 2/5 ∞ R 3 f (t, x, v)|v| 2 dv 3/5
and after integration we obtain the claim.

B Uniform lower bound on n e

The purpose of this section is to prove several properties that we have used throughout the paper for the elliptic equation

     -λ 2 ∆ ln n e + n e = n I , x ∈ Ω, ∂n e ∂ν = 0, x ∈ ∂Ω,
with a right hand side data satisfying n I ≥ 0, n I = M initial > 0 and n I ∈ L 5/3 .

We are going to prove the estimate

Lemma 2. Let n I ∈ L ∞ (0, T ; L 5/3 ∩ L 1 (Ω))
and n e a strong solution to the above equation, then we have the two-sided control 0

< K -( n I 5/3 ) ≤ n e ≤ K + ( n I 5/3 ), (35)
for some continuous positive functions K ± (•) with K + > 1 increasing, K -decreasing.

We also recall that integration of the equation gives the electric neutrality relation

n e 1 = n I 1 .
We can now explain why the total energy ( 8) is well defined for weak solutions and also at initial time. From the lower and upper bound in (35), we conclude that n e ln n e ∈ L 1 (Ω). Also, multiplying the equation by ln n e , we find λ

Ω |∇ x ln n e | 2 dx = Ω (n I -n e ) ln n e ≤ 2 ln K + n I 1 .
Proof. Then, we argue in two steps. Firstly, we multiply the equation by n Secondly we use the elliptic regularity theory to conclude that ln n e -ln n e Ω ∈ W 2,p , 1 ≤ p ≤ 5/3 and thus

ln n e -ln n e Ω ∈ L q , ∀q > 1, 1 q = 1 p - 2 3 ,
where φ Ω denotes the average of the L 1 (Ω) function φ over Ω. Finally, because 2 3 > 3 5 we conclude from the Morrey estimates [START_REF] Evans | Partial Differential Equations[END_REF] that ln n e -ln n e Ω ∞ ≤ C( n I 5/3 ).

The result follows immediately thanks to the control of ln n e Ω through (35).

C Compactness of the magnetic field

We have also used the Sobolev injection for Maxwell equations and we recall it in this appendix. We introduce the following spaces

H curl (Ω) = {b ∈ L 2 (Ω) 3 / ∇ ∧ b ∈ L 2 (Ω) 3 }, H div (Ω) = {b ∈ L 2 (Ω) 3 / ∇ • b ∈ L 2 (Ω) 3 }, X N (Ω) = {b ∈ H curl (Ω) ∩ H div (Ω) / b ∧ n = 0 on ∂Ω}. (36)
We recall the following result (see [START_REF] Dautray | Mathematical analysis and numerical methods for sciences and technology[END_REF])

Theorem 2. Assume that the domain Ω is of class C 1,1 . Then the space X N (Ω) is continuously imbedded in H 1 (Ω) 3 .

D Kinetic averaging lemma

We recall here one result of the theory of averaging lemmas for kinetic equations. When f (t, x, v) is solution of a kinetic equation, it cannot be more regular that the initial data or the right hand-side. However, averages in velocity gain regularity. Recall that the macroscopic quantity f ψ is defined as

f ψ (t, x) = R d f (t, x, v)ψ(v)dv
where ψ is a given function in C ∞ c (R d ) (i.e. smooth with compact support), the averaging lemmas aim at proving compactness properties on f ψ . The first version of these averaging Lemma has been established by Golse, Lions Perthame and Sentis [START_REF] Golse | Regularity of the moments of the solution of a transport equation[END_REF] for f ∈ L 2 t,x,v solution of the equation ∂ t f + v • ∇ x f = S with S ∈ L 2 t,x,v : it was proved that locally in time f ψ ∈ H 1/2 t,x . This version has been then be extended (by complex interpolation) for the L p framework, 1 < p < ∞ in [START_REF] Diperna | Lp regularity of velocity averages[END_REF]. Moreover, more complex versions have been proved by Di Perna, Lions and Meyer [START_REF] Diperna | Lp regularity of velocity averages[END_REF], who treat the case where S is the k-th derivative in velocity with a fractional derivative in x strictly less than one.

We use here a version proved by Perthame and Sougadinis [START_REF] Perthame | A limiting case for velocity averaging[END_REF] (the equality is the exponent is due to Bouchut [START_REF] Bouchut | Kinetic equations and asymptotic theory[END_REF]) which expresses an optimal gain of regularity (a full derivative) Theorem 3. Let 1 < q < ∞ and f, g = (g 1 , ..., g d ) belong to L q t,x,v (R 1+3+3 ) and satisfy

∂f ∂t + v • ∇ x f = ∂ ∂v . (I -∆ t,x,v ) 1/2 g . ( 37
)
Then f ψ ∈ L q t,x (R 1+3 ) and there exists C(p, ψ) such that f ψ q ≤ C(q, ψ) f 1-α q g α q for a positive exponent α ≤ 1 2 min( 1 q , 1 -1 q ). In fact this theorem also proves strong compactness, which is the way we use it in our context. Consider a truncation function χ = χ(t, x) with suppχ ⊂ (0, T ] × Ω which is fixed in this paragraph. Define the sequence of functions w ε as w ε = χF ε f ε and z ε = Zf ε , with Z = ∂χ ∂t + v.∇ x χ. They satisfy in R 1+3+3 ∂ ∂t

(χf ε ) + v • ∇ x (χf ε ) = ∂ ∂v .w ε + z ε . (38)
We know that w ε is bounded in L 1 ∩ L p t,x,v (R 1+3+3 ) for some p > 1 (here p = 30/29) and z ε and χf ε are bounded in L 1 ∩ L ∞ t,x,v (R 1+3+3 ). We are going to prove the following result for an exponent q (with 1 < q < p). Lemma 3. Consider, after extraction, the weak limit f * in L p t,x,v (R 1+3+3 ) of the sequence f ε weakly then

ρ ε ψ = χ f ε ψ → ρ * ψ = χ f * ψ , strongly in L q t,x,v (R 1+3+3 ). (39)
Proof. We first define g ε = (I -∆ x,t,v ) -1/2 w. ε

By regularizing effects the family g ε is compact in L 1 t,x,v (R 1+3+3 ) : indeed g ε is given by a convolution product between w ε and the fundamental solution of (I -∆ x,t,v ) -1/2 . Moreover g ε is bounded in L p t,x,v (R 1+3+3 ) and then, by interpolation argument, it is compact in L q t,x,v (R 1+3+3 ) for some q (with 1 < q < p). Since the family {χf ε } is uniformly bounded in L q t,x,v as well as z ε , the theorem recalled above (the term z ε can be written as a v -derivative without loss of generality) applied to χf ε allows to claim that the following bound holds ρ ε ψ q ≤ C(q, ψ)C g ε α q .

Then, there exists a subsequence f ε , w ε and functions f * , ρ * ψ , g * , w * such that

           f ε ⇀ f * , w ε ⇀ w * weakly in L p t,x,v (R 1+3+3 ), ρ ε ψ ⇀ ρ * ψ
weakly in L q t,x,v (R 1+3+3 ), g ε → g * strongly in L q t,x,v (R 1+3+3 ).

and g * = (I -∆ x,t,v ) -1/2 w * . Of course we have also, passing to the weak limit,

∂ ∂t + v • ∇ x )(χf * = ∂ ∂v • w * + Zf * = ∂ ∂v • ((I -∆ x,t,v ) 1/2 g * ) + Zf * .
We combine this with (38) and we get

( ∂ ∂t + v • ∇ x )[χf ε -χf * ] = ∂ ∂v • ((I -∆ x,t,v ) 1/2 (g ε -g * )) + Z(f ε -f * ).
And according to the previous theorem (the term Z(f ε -f * ) can be written as a v-derivative also) once again we see that χ f ε ψ -χ f * ψ q ≤ C(q, ψ) • g εg * α q f ε -f * (1-α) q .

That is to say, the property (39) holds.

  ∂f ∂t + v.∇f + ∂ ∂v .[(-T e ∇ ln n e +v ∧ B) f ] = 0, -λ 2 ∆ln n e = f (v)dv-n e ,

4 5 ≤ 1 + y 13 10

 4113 for positive y, we get the more compact form R ≤ δ 1 + δ 2 E pert tot 13 10

n

  I u I = n I u I (t k ) + t t k ∇ ∧ B(s)ds ∧ d. It shows that n I u I integro-differential and linear with respect to B. Plugging this form of n I u I in equation b), we end up with a linear equation for B. This linear integro-differential equation is well posed under general assumptions.

+ n 5/ 3 e= n I n 2/ 3 e≤ n I 5 /3 n e 2 / 3 5/ 3 ,

 335233 from which we conclude the bound n e 5/3 ≤ n I 5/3 .


Fréderique Charles

Cédrick Copol

Stéphane Dellacherie

Jean-Marc Mounsamy

Numerical simulation by a random particle method of Deuterium-Tritium fusion reactions in a plasma

We propose and we justify a Monte-Carlo algorithm which solves a spatially homogeneous kinetic equation of Boltzmann type that models the fusion reaction between a deuterium ion and a tritium ion, and giving an α particle and a neutron. The proposed algorithm is validated with the use of explicit solutions of the kinetic model obtained by replacing the fusion cross-section by a Maxwellian cross section.

Résumé

On propose et on justifie un algorithme de type Monte-Carlo permettant de résoudre un modèle cinétique homogène en espace de type Boltzmann modélisant la réaction de fusion entre un ion deutérium et un ion tritium, et donnant une particule α et un neutron. L'algorithme proposé est par ailleurs validé via des solutions explicites du modèle cinétique obtenues en remplaçant la section efficace de fusion par une section efficace maxwellienne.

Introduction

One of the goals of the futur International Thermonuclear Experimental Reactor is to obtain fusion reactions between two isotopes of hydrogen in a confined plasma. We study here one of the possible reaction: T(d,n) 4 He, between Deuterium (D) and Tritium (T ) ions producing an α particle and a neutron (n) with the following scheme:

D + T → α + n + 17, 59 MeV.
When two particles of Deuterium and Tritium fuse, they disappear and produce an α particule and a neutron with a mass default ∆M . According to Einstein's law, the energy produced by this mass 1 is Q := ∆M c 2 = 17, 59 MeV (where c is the light speed). Here we assume that this energy is fully converted into kinetic energy. The very energetical α particles produced have a significant role in ITER: they take part in the heating of the plasma, contribute to the evolution of the electrical field and can cause instabilities. For these reasons, the study of their trajectories and their interactions with other species is an important point for the comprehension of instabilities in a fusion plasma. In [START_REF] Farengo | Simulation of alpha particle current drive and heating in spherical tokamaks[END_REF] for example, the α particles current drive and heating is studied, using an expression of the fusion reaction rate τ which is the rate of α particles produced per unit volume per second (see §1.2). However, in areas where the plasma may not be at thermodynamical equilibrium, it is more appropriated -in particular because classical expressions of the nuclear fusion rate are not true -to use a kinetic modelling and to compute the velocities distribution of the α particles. In [START_REF] Dellacherie | Nuclear Collisions Models with Boltzmann Operators[END_REF], Dellacherie and Sentis introduced a kinetic model to describe nuclear collisions. In a spatially homogeneous context, the model writes:

                             ∂f D ∂t = -P D (f D , f T ), (a) ∂f T ∂t = -P T (f D , f T ), (b) ∂f α ∂t = G α (f D , f T ), (c)
∂f n ∂t = G n (f D , f T ), (d) (1)
where f a : R + × R 3 → R + is a non negative function, for a ∈ {D, T, α, n} ; f a (t, v a ) corresponds to the number densities1 in ions Deuterium, ions Tritium, neutrons, and α particles respectively at time t ∈ R + at velocity v a ∈ R 3 . In [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF], operators P a (f D , f T ), a ∈ {D, T } are loss operators, which modelize the disappearing of a D or T particle after a collision, wheras G a (f D , f T ), a ∈ {α, n} are gain operators, corresponding to the creation of a neutron or an α particle. An asymptotic analysis of the gain operator when the energy of the reaction Q → +∞ is also performed, and leads to an approched model where the gain operators are substituted by Dirac masses on S 2 . We propose here a random particle method simulation of the model [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF], and we compare this last model with another model obtained as a limit when the ratio between the kinetic temperature of the plasma and Q tends to 0. The paper is organized as follows. In the first Section, we present the collision model introduced by Dellacherie and Sentis in [START_REF] Dellacherie | Nuclear Collisions Models with Boltzmann Operators[END_REF] and the asymptotic model. We also establish expression of explicite solutions in the case of a Maxwellian cross-section of collision. In the second Section, we present the numerical method we developped to solve the model (1). This method is based on a Monte-Carlo simulation of the loss and Gain operators, and is close from the Nanbu's method for Boltzman operator. The third Section is devoted to numerical results. First, we validate our algorithm thanks to explicit solutions obtained in the case of the Maxwellian cross-section. Then, we present the numerical simulation of model ( 1) with the fusion cross-section. Finally, we investigate the behaviour of the function density in energy of α particles when the kinetic temperature of the plasma decreases, in order to discuss the oportunity of using an asymptotic model.

Kinetic modeling of D-T fusion reaction 1.Collision operators

The model established in [START_REF] Dellacherie | Nuclear Collisions Models with Boltzmann Operators[END_REF] is based on the following hypothesis: -the energy corresponding to the mass default Q = ∆M c 2 is totally converted into kinetic energy on α particles and neutron, -we neglect the mass default ∆M with respect to the masses of particles which means that we assume that

m D + m T = m α + m n (2)
where m D , m T , m α and m n are the masses of ions Deuterium, Tritium, α and of a neutron respectively, -the momentum is conserved during a collision. Consequently, the momentum conservation and the energy balances write

m D v D + m T v T = m α v α + m n v n , (a)
m D v 2 D + m T v 2 T + 2Q = m α v 2 α + m n v 2 n . (b) (3)
Let us define the following notations:

µ DT = m D m T m D + m T , µ α n = m α m n m α + m n , ζ = 2Q/µ αn . (4)
From (3), we deduce that there exists Ω ∈ S2 such that the velocities of the created particles v α and v n write

v α = v α (v D , v T , Ω) = m D v D + m T v T m D + m T + Ω µ αn m α ζ 2 + µ DT µ αn |v D -v T | 2 , (a)
v n = v α (v D , v T , Ω) = m D v D + m T v T m D + m T -Ω µ αn m n ζ 2 + µ DT µ αn |v D -v T | 2 . (b) (5)
It is then possible to derive the expression of the loss and gain operators. The loss operator P D writes (we have a similar expression for P T ):

P D (f D , f T )(v D ) = f D (v D ) ˆR3 ˆS2 f T (v T )q(|v D -v T |, Ω)dv T dΩ, (6)
and the gain operator G α writes on a weak form 2 (we have a similar expression for G n ):

ˆR3 G α (f D , f T )(v α )ϕ(v α )dv D = ˆR3 ˆR3 ˆS2 f D (v D )f T (v T )ϕ(v α )q(|v D -v T |, Ω)dv D dv T dΩ, (7)
where v α depends on v D , v T and Ω through [START_REF] Farengo | Simulation of alpha particle current drive and heating in spherical tokamaks[END_REF]. The collision kernel q is here isotrope -that is

q(|v D -v T |, Ω) = q(|v D -v T |)
-and is defined by

q(|v D -v T |) = |v D -v T |σ(|v D -v T |),
where σ(|v Dv T |) is called the collision cross-section. We will here consider two types of collision cross-section:

-the nuclear cross-section which has the following form3

σ f (|v D -v T |) = p 5 + p 2 /[(p 4 -p 3 E) 2 + 1] E[exp(p 1 / √ E) -1] (8)
where [START_REF] Duane | Fusion cross-section theory[END_REF][START_REF] Huba | NRL plasma formulary[END_REF], -the Maxwellian cross-section 4 given by

E = 1 2 m D |v D -v T | 2 and p i ∈ R + are the Duane coefficients
σ m (|v D -v T |) = λ 4π|v D -v T | ( 9
)
where λ is a constant (in m 3 •s -1 ). Thus, we have q(|v Dv T |) = λ/4π. We use this cross-section to obtain explicit solutions which enables us to validate the numerical method. Remark 1. It should be possible to prove the existence and uniqueneess of non-negative solutions of (1) with initial conditions with appropriated assumptions on q and on the initial condition in the same spirit as in [START_REF] Charles | Modélisation mathématiques et étude numérique d'un aérosol dans un gaz raréfié. Application à la simulation du transport de particules de poussière en cas d'accident de pert de vide dans ITER[END_REF]. However, we do not study this problem here and we assume that the moments of order 0, 1 and 2 in velocity of the density functions f a (a ∈ {D, T, α, n}) are well defined.

Figure 1: Cross-section σ f (in Barns, 1 Barn=10 -28 m 2 ) according to E (in keV) in log-log scale.

Definition of some physical quantities

We define here some physical quantities used in the following. For each species a ∈ {D, T, α, n}, we define the number density ρ a (without unit), the macroscopic velocity u a (in m•s -1 ) and the kinetic temperature θ a (here in Joules 5 ) with

           ρ a (t) = f a (t, •) , u a (t) = 1 ρ a f a (t, •)v if ρ a (t) > 0, θ a (t) = m a 3ρ a f a (t, •)(v -u a ) 2 if ρ a (t) > 0, (10)
where g represents the average in velocity of a function g ∈ L 1 v , that is ´R3 g(v)dv. We also define the fusion reaction rate which corresponds to the number of particles α (or of neutrons) created per unit of volum and of time:

τ (t) = τ (f D , f T )(t) := ˆR3 ˆR3 ˆS2 f D (t, v D )f T (t, v T )q(|v D -v T |, Ω)dΩdv T dv D , (11)
the total momentum W :

W (t) = m D ρ D u D + m T ρ T u T + m α ρ α u α + m n ρ n u n , (12)
and the total kinetic energy E tot :

2E tot (t) = m D f D (t, v D )v 2 D + m T f T (t, v T )v 2 T + m α f α (t, v α )v 2 α + m n f n (t, v n )v 2 n . (13)
Let us underline that τ corresponds to

P a (f D , f T ) (a ∈ {D, T }) and to G a (f D , f T ) (a ∈ {α, n}).
Moreover if f a , a ∈ {D, T, α, n} are solutions of (1) then, thanks to relation (3), we have:

             τ = - d dt ρ D = - d dt ρ T = d dt ρ α = d dt ρ n , (a)
dW dt = 0, (b)
dE tot dt = Qτ (t). (c) (14)
5. We can also define in Kelvin the temperature Ta of the particles of type a by using the relation θa = k B Ta where k B = 1.3806 • 10 -23 J•K -1 is the Boltzmann constant.

An asymptotic model

Let us introduce the thermal velocities V th a := 2θ a /m a (a ∈ {D, T }), V th p = max(V th D , V th T ) and

ε = V th p ζ . (15)
For θ D = θ T = 10 keV 6 (which corresponds to a temperature in K of T D = θ D /k B = 1.16 • 10 8 K, which is the classical order of magnitude of the temperature in a Tokamak plasma), we have ε = 0.015. In this situation, we see in (5) that, formally, we have

|v α | ≈ V α := µ αn m α ζ and |v n | ≈ V n := µ αn m n ζ.
Then, the velocities of created particles α and n are very close to the standard creation velocities V α and V n , which correspond to the well-know values of energy

     E α = 1 2 m α V α 2 = 3.52 MeV, E n = 1 2 m n V n 2 = 14.08 MeV. ( 16
) More precisely, one can prove that if q ∈ L ∞ (R), f D and f T are in L 1 (R 3 ; (1 + |v| 2 )dv), the term G a (f D , f T ) G a (f D , f T ) (a ∈ {α, n})
converges weakly in the sens of measures to the Dirac mass on

S 2 1 4π δ a defined by ˆR3 ϕ(v)δ a (v)dv = ˆS2 ϕ(V a Ω)dΩ, ∀ϕ ∈ C(R 3 ). ( 17
)
When ε ≪ 1 in (1), we can therefore consider an approached model by replacing the gain operator G a by its weak limit, which gives

                         ∂f D ∂t = -P D (f D , f T ), (a) ∂f T ∂t = -P T (f D , f T ), (b)
∂f α ∂t = 1 4π τ (f D , f T ) δ α , (c)
∂f n ∂t = 1 4π τ (f D , f T ) δ n , (d) (18)
where τ (f D , f T ) is defined by [START_REF] Li | Fusion cross-sections for inertial fusion energy[END_REF].

Remark 2. On the contrary of model (1), density functions solutions of model (18) does not verify equations (14)(b,c).

Explicit solutions in the case of a Maxwellian cross-section

Proposition 1. We consider the system (1) with the following initial conditions

f D (0, v) = f ini D (v) ≥ 0, f T (0, v) = f ini T (v) ≥ 0, f α (0, v) = 0, f n (0, v) = 0, ( 19
)
where the cross-section σ in the operators P D , P T , G α and G n is given by the Maxwellian cross-section σ m (9). When 0 < ρ D (0) = ρ T (0) =: ρ ini (where ρ D and ρ T are defined by (10)), f D and f T are given by

         f D (t, v D ) = f ini D (v D ) λρ ini t + 1 , f T (t, v T ) = f ini T (v T ) λρ ini t + 1 (20) 6. 1 keV=1.602 • 10 -16 J
and f α and f n are given by the weak forms

f α (t, •)ϕ = λ(ρ ini ) 2 t λρ ini t + 1 R α (f ini D , f ini T )(ϕ), f n (t, •)ϕ = λ(ρ ini ) 2 t λρ ini t + 1 R n (f ini D , f ini T )(ϕ), (21)
for all ϕ such that f a ϕ (a ∈ {α, n}) is well defined, R a (a ∈ {α, n})) being given by

R a (f ini D , f ini T )(ϕ) = 1 4π(ρ ini ) 2 ˆR3 ˆR3 ˆS2 f ini D (v D )f ini T (v T )ϕ(v a )dv D dv T dΩ, a ∈ {α, n}, (22)
where

v a = v a (v D , v T , Ω
) is given by (5).

We deduce from Proposition 1 that the fusion reaction rate τ is given by

τ (t) = λ(ρ ini ) 2 (λρ ini t + 1) 2 , ( 23
)
that the number densities are given by

ρ D (t) = ρ T (t) = ρ ini λρ ini t + 1 , ρ α (t) = ρ n (t) = λ(ρ ini ) 2 t λρ ini t + 1 for t ≥ 0, (24)
that the macroscopic velocities are given by

u D (t) = u D (0), u T (t) = u T (0), u a (t) = R a (f ini D , f ini T )(1 → v) for t > 0, a ∈ {α, n}, (25)
and that the kinetic temperatures are given by

θ D (t) = θ D (0), θ T (t) = θ T (0), θ a (t) = m a 3 R a (f ini D , f ini T )(1 → (v -u a ) 2 ) for t > 0, a ∈ {α, n}. (26
) Let us remark that the time (λρ ini ) -1 defines an order of magnitude of the scale of time of the fusion reaction. Moreover, we have lim

t→+∞ ρ D (t) = lim t→+∞ ρ D (t) = 0 and lim t→+∞ ρ α (t) = lim t→+∞ ρ n (t) = ρ ini .
We can also obtain explicit solution in the case of model (18). Proposition 2. We consider system (18) with initial conditions (19) where the cross-section σ in the operators P D , P T is given by the Maxwellian cross-section σ m [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]. Then, f α and f n are given by

f a (t, v a ) = λ(ρ ini ) 2 t λρ ini t + 1 • δ a (v a ) 4π , a ∈ {α, n}, (27)
and the macroscopic quantities are given by:

ρ α (t) = ρ n (t) = λ(ρ ini ) 2 t λρ ini t + 1 for t ≥ 0, ( 28
)
u α (t) = u n (t) = 0, θ α (t) = m α 3 V α 2 , θ n (t) = m n 3 V n 2 for t > 0. ( 29
)

Numerical method

The numerical resolution of systems (1) and ( 18) is based on a particle method coupled to a Monte-Carlo simulation of fusion operators, like in Nanbu's method [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF][START_REF] Pareschi | An introduction to Monte Carlo methods for the Boltzmann equation[END_REF][START_REF] Perthame | Introduction to the theory of random particle methods for Boltzmann equation[END_REF]. The initial densities f ini D and f ini T are approached by the linear combination of Dirac masses

f ini a (v) ≈ f 0 a (v) := N 0 a k=1 ω a,k δ(v -V 0 a,k ) for a ∈ {D, T }, (30)
which corresponds to represent the distribution of particles of species a by N 0 a numerical particles having a numerical weight ω a,k and a velocity V 0 a,k . In order to use a symetrical method, we choose the same numerical weight ω a,k for all numerical particles and for both species D and T . The number N 0 D being fixed, we then take

ω = ρ ini D N 0 D and N 0 T = ρ ini T ω . (31)
The velocities V 0 a,k of numerical particles are computed thanks to a Monte-Carlo method following the initial density of probability f ini a /ρ ini a . The derivative in time is approched by an explicit Euler scheme, and the computation of gain and loss operators is based on a Monte-Carlo simulation. Moreover, the initial densities f ini α and f ini n are supposed to be equal to zero. Finally, the densities f n a (a ∈ {D, T, α, n}) at any time t n are given by

f n a (v) = ω N n a i=1 δ(v -V n a,i ) for a ∈ {D, T, α, n} (32)
where N n a and (V n a,i ) {1≤i≤N n a } are deduced from the Monte-Carlo simulation. In the sequel, we detail the Monte-Carlo algorithm, firstly, for any cross-section such that q ∈ L ∞ (R + ) (the fusion cross-section given by ( 8) -see also Figure 1 -satifies q ∈ L ∞ (R + )) and, secondly, for a Maxwellian cross-section [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF]. We only present the algorithm for the densities f D and f α , the principle being the same for f T and f n .

Fusion cross-section

2.1.1 Other expressions of P D (f D , f T ) and G α (f D , f T ) The computation of the terms P a (f n D , f n T )(v) and G a (f n D , f n T )(v)
with a fusion cross-section is based on the following propositions (see Annex A for the proofs): Proposition 3. Let ϕ be a test function (for example ϕ ∈ C c (R 3 )) and let P D be the loss operator given by [START_REF] Goncharov | Practical Calculation of Nuclear Fusion Power for a Toroidal Plasma Device with Magnetic Confinement[END_REF] where the collision kernel q is in L ∞ (R + ). Let q max be an upper bound of q. We set

f n a (v) = ω N n a i=1 δ(v -V n a,i ) for a ∈ {D, T }.
We assume that

∆t := t n+1 -t n ≤ 1 λ max ω min 1 N n T , 1 N n D , with λ max = 4πq max . (33)
Then, we have

∆t ˆR3 ϕ(v)P D (f n D , f n T )(v)dv = η n ω ♯P P N n k=1 ϕ(V n D,r k ) ˆ1 0 1 {0≤u≤p n r k ,s k } du (34)
where

             N n = ⌊N n D N n T ω∆tλ max ⌋, η n = N n D N n T ω∆tλ max N n , p n r k ,s k = q(|V n D,r k -V n T,s k |) q max , ( 35
)
and P is the set of possible ways to make

N n distinct pairs (r k , s k ) {1≤k≤N n } among {1, . . . , N n D } × {1, . . . , N n T }.
Proposition 4. Let G α be the gain operator defined by [START_REF] Heltemes | Analysis of an Improved Fusion Reaction Rate Model for Use in Fusion Plasma Simulations[END_REF]. Under assumptions of Proposition 3, we have

∆t ˆR3 ϕ(v)G α (f n D , f n T )(v)dv = η n ω ♯P P N n k=1 ˆS2 ˆ1 0 ϕ(V n α,r k ,s k ,Ω )1 {0≤u≤p n r k ,s k } du dΩ 4π ( 36
)
where

V n α,r k ,s k ,Ω = v α (V n D,r k , V n T,s k , Ω
) is given by formula (5)(a).

Let us note that, typically, N n is of the order of 10 3 . Thus, we can approximate η n with 1.

Monte-Carlo approximation

The terms that we have obtained are then computed by a Monte-Carlo procedure. More precisely, we approach

ω ♯P P N n k=1 ϕ(V n D,r k ) ˆ1 0 1 {0≤u≤p n r k ,s k } du by ω N n k=1 ϕ(V n D,r k )1 {0≤u k ≤p n r k ,s k } ( 37
)
and

ω ♯P P N n k=1 ˆS2 ˆ1 0 ϕ(V n α,r k ,s k ,Ω )1 {0≤u≤p n r k ,s k } du dΩ 4π by ω N n k=1 ϕ(V n α,r k ,s k ,Ω k )1 {0≤u k ≤p n r k ,s k } (38)
where the set of pairs (r k , s k ) {1≤k≤N n } is chosen randomly among P with an uniform law, u k is randomly computed with an uniform law on [0, 1] and Ω k is randomly computed with an uniform law on S 2 . Let us denote P D and G α the approximated operators P D and G α that we obtain by using respectively approximations (37) and ( 38), which means that (by approximating η n with 1)

               ∆t ˆR3 ϕ(v) P D (f n D , f n T )(v)dv = ω N n k=1 ϕ(V n D,r k )1 {0≤u k ≤p n r k ,s k } , (a)
∆t

ˆR3 ϕ(v) G α (f n D , f n T )(v)dv = ω N n k=1 ϕ(V n α,r k ,s k ,Ω k )1 {0≤u k ≤p n r k ,s k } . (b) (39)
Moreover, we use an explicit Euler scheme for the time discretization of equations [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF]. Then, if f n D and f n α are the approximation at time t n of the exact densities f D (t n , •) and f α (t n , •), we define

f n+1 D and f n+1 α by    f n+1 D = f n D -∆t P D (f n D , f n T ), (a) f n+1 α = f n α + ∆t G α (f n D , f n T ). (b) (40)
Thus, by using (32), (39)(a) and (40)(a), we obtain

ˆR3 ϕ(v) f n+1 D (v)dv = ω N n D i=1 ϕ(V n D,i ) -ω N n k=1 ϕ(V n D,r k )1 {0≤u k ≤p n r k ,s k } = ω N n D i=1 ϕ(V n D,i ) -ω k∈A ϕ(V n D,r k ) = ω k∈{1,...,N n D }\A ϕ(V n D,r k )
where we set

A = {k ∈ {1, . . . , N n }/ 0 ≤ u k ≤ p n r k ,s k } and N n+1 D = ♯({1, . . . , N n D } \ A). (41)
Thus, we have

ˆR3 ϕ(v) f n+1 D (v)dv = ω N n+1 D i=1 ϕ(V n+1 D,i ) that is to say f n+1 D (v) = ω N n+1 a i=1 δ(v -V n+1 a,i )
where

{V n+1 D,1 , . . . , V n+1 D,N n+1 D
} is the set of velocities obtained when we remove

{V n D,r1 , . . . , V n D,r N n } from {V n D,1 , . . . , V n D,N n D }.
In the same way, we deduce from (32), (39)(b) and (40

)(b) that ˆR3 ϕ(v) f n+1 α (v)dv = ω N n α i=1 ϕ(V n α ) + ω N n k=1 ϕ(V n α,r k ,s k ,Ω k )1 {0≤u k ≤p n r k ,s k } = ω N n+1 α i=1 ϕ(V n+1 α,i )
that is to say

f n+1 α (v) = ω N n+1 α k=1 δ(v -V n+1 α,k )
where

N n+1 α = N n α + ♯A and    V n+1 α,i = V n α,i if 1 ≤ i ≤ N n α , {V n+1 α,i , N n α < i ≤ N n+1 α } = {V n α,r k ,s k ,Ω k , k ∈ A}.

Algorithm

We therefore deduce the following algorithm:

Algorithm 1. Random Particle method for the model (1) with a fusion cross-section.

At time step t n , we perform the following operations:

i) We compute N n = ⌊N n D N n T ω∆tλ max ⌋. The times step ∆t = t n+1 -t n (which could depend of n) is such that condition (33) is satisfied. Thus, N n is smaller than N n D and N n T . If N n ≪ min(N n D , N n T ), one can increase ∆t.
ii) We select successively N n distinct numerical particles of D (respectively of T ) the k th particle being randomly chosen with an uniform discret law among the N n D -k +1 (respectively N n T -k +1) remaining particles. This corresponds to draw the pairs (r k , s k ) k∈{1,...,N n } with an uniform law among the ♯P possible set of N n pairs. iii) For each of the N n pairs of velocities (V D,r k , V T,s k ), we draw u k following an uniform law on [0, 1] and we compute the term

p r k ,s k = q(V D,r k -V T,s k ) q max ,
where q(ξ) = |ξ|σ f (|ξ|), σ f being given by [START_REF] Huba | NRL plasma formulary[END_REF]. If p r k ,s k > u k , the two particles of velocities V D,r k and V T,s k interact (in other word, we have k ∈ A, where A is defined by (41)). iv) For each pair of interacting particles (that is to say for each k ∈ A), we draw 7 Ω k with an uniform law on S 2 , and V α,r k ,s k ,Ω k and V n,r k ,s k ,Ω k are computed thanks to formulae [START_REF] Farengo | Simulation of alpha particle current drive and heating in spherical tokamaks[END_REF].

v) The ♯A pairs of interacting numerical particles of D and T are removed, and ♯A pairs of particles of α and n particles are created with the previous velocities

V α,r k ,s k ,Ω k and V n,r k ,s k ,Ω k , 1 ≤ k ≤ N n .
The velocities of other particles remain unchanged.

Remark 3. Since a numerical fusion reaction represents ω real fusion reaction, the total number of fusion reactions during a time step ∆t is ω

N n k=1 1 {0≤u k ≤pr k ,s k } which corresponds to a Monte-Carlo approximation of ω ♯P P N n k=1 ˆ1 0 1 {0≤u≤pr k ,s k } du = N n N n D N n T N n D i=1 N n T i=1 ω q(V n D,i , V n T,j ) q max ≃ η n N n N n D N n T N n D i=1 N n T i=1 ω q(V n D,i , V n T,j ) q max (since η n ≃ 1) = 4π∆t ˆR3 ˆR3 f n D (v D )f n T (v T )q([v D -v T |)dv D dv T = ∆t τ (f n D , f n T ) ≃ ∆ tτ (t n ),
where τ (t n ) is the exact fusion reaction rate given by (11) at time t n .

Remark 4. Algorithm 1 preserves exactly the total number of particles and the total momentum W . Moreover, the total energy verifies

E n+1 tot = E n tot + Q∆t τ (f n D , f n T ).
In other words, ( 14) is satisfied at the discrete level.

Maxwellian cross-section

The numerical method introduced in §2.1 for any cross-section can be simplified for a Maxwellian cross-section. Indeed, by taking q max = q = λ 4π , we have p i,j = 1 for all couple of velocity (V D,i , V T,j ), and then 1 {0≤u k ≤pr k ,s k } = 1 for all k ∈ {1, . . . , N n }. Consequently, there is no need to draw u k . Thus, denoting f n D and f n α the approximation of the exact density f D (t n , •) and f α (t n , •) at time t n , we approach the exact densities at time t n+1 by f n+1 D and f n+1 α with

f n+1 D (v) = ω N n+1 D k=1 δ(v -V n+1 D,k ), f n+1 α (v) = ω N n+1 α k=1 δ(v -V n+1 α,k ) (42)
where

{V n+1 D,1 , . . . , V n+1 D,N n+1 D } is the set of velocities obtained when we remove {V n D,r1 , . . . , V n D,r N n } from {V n D,1 , . . . , V n D,N n D }, N n+1 α = N n α + N n and    V n+1 α,k = V n α,k for 1 ≤ k ≤ N n α V n+1 α,N n α +k = V n α,r k ,s k ,Ω k for 1 ≤ k ≤ N n .
This leads to the following simplified algorithm:

Algorithm 2. Random Particle method for the model (1) with Maxwellian cross-section:

At time step t n , the steps (i) and (ii) are the same as in Algorithm 1. Steps (iii)-(v) are replaced by:

7. We can set Ω k = (a k cos(2πǫ k ), a k sin(2πǫ k ), b k ) with b k = 1 -2u k , a k = 1 -b 2 k , u k and ǫ k being computed following an uniform law on [0, 1].
iv) For each of the N n pairs of velocities (V D,r k , V T,s k ) selected, we compute Ω k with an uniform law on S 2 , and V α,r k ,s k ,Ω k and V n,r k ,s k ,Ω k are computed thanks to formulae [START_REF] Farengo | Simulation of alpha particle current drive and heating in spherical tokamaks[END_REF].

v) The N n selected pairs of numerical particles of D and T are removed, and N n pairs of α and n particles are created with the previous velocities

V α,r k ,s k ,Ω k and V n,r k ,s k ,Ω k , 1 ≤ k ≤ N n .
The velocities of other particles remain the same.

Asymptotic model

We now consider the system (18) with a collision cross-section q ∈ L ∞ in the operators P a (f D , d T ) (a ∈ {D, T }). This system is an approximation of system [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF]. The numerical resolution of equations (18)(a,b) does not change since (18)(a,b) and ( 1)(a,b) are identical. For the resolution of (18)(c,d), we write thanks to definition (17) of δ a the equation on f a (a ∈ {α, n}) on the following weak form:

d dt f a (t, •)ϕ = τ (f D , f T )(t) ˆS2 ϕ(V a Ω) dΩ 4π . ( 43
)
Similar computations to those of Proposition 4 show that under condition (33) on the time step ∆t, we can approach by a Monte-Carlo procedure

∆t τ (f n D , f n T ) ˆS2 ϕ(V a Ω) dΩ 4π by η n ω N n k=1 ϕ(V a Ω k )1 {0≤u k ≤pr k ,s k }
where u k , r k , s k are the random variables drawn for the numerical simulation of equation ( 1)(a,b) (steps (ii) and (iii) of Algorithm 1), and where (Ω k ) k∈{1,...,N n } are drawn following an uniform law on S 2 . Thus, by applying as previously an Euler scheme for the time discretization, we obtain (by using the fact that η n ≃ 1)

f n+1 a ϕ = f n a ϕ + ω N n k=1 ϕ(V a Ω k )1 {0≤u k ≤pr k ,s k } . (44)
Finally, we obtain

f n+1 a (v) = ω N n+1 a k=1 δ(v -V n+1 a,k ) for a ∈ {α, n} with N n+1 a = N n a + ♯A and    V n+1 a,i = V n a,i for 1 ≤ i ≤ N n a , {V n+1 a,i , N n a < i ≤ N n+1 a } = {V a Ω k , k ∈ A},
A being defined by (41). This leads to the following simplified algorithm: Algorithm 3. Random Particle method for the model (18) with fusion cross-section.

At time step t n , steps (i) and (iii) are the same as in Algorithm 1. Steps (iv)-(v) are replaced by: iv) The ♯A pairs of interacting numerical particles of D and T are removed. The velocities of other particles of species D and T remain the same.

v) We draw ♯A vectors of S 2 Ω 1 , . . . Ω ♯A uniformly on S 2 , and ♯A pairs of particles of α and n are created with velocities {V a Ω k , k ∈ A}. The velocities of other particles α and n remain the same.

Remark 5. We can easily adapt Algorithm 3 for a Maxwellian cross-section: the step (iii) is removed, and we replace ♯A by N n in steps (iv) and (v).

Remark 6. We see that the model (18) is as costely to compute as model (1), and therefore does not present any numerical interest. However, if we are only interested in the evolution of f α and f n, we can introduce the following model

     ∂f α ∂t = 1 4π τ (t) δ α , ∂f n ∂t = 1 4π τ (t) δ n (45)
where τ (t) is a given fusion reaction rate. Then, at each time step, we only need to compute the step (v) of Algorithm 3 (with ⌊ ∆t τ (t n ) ω ⌋, ω being the numerical weight we choose for numerical particles, instead of ♯A). The numerical cost is lower than the one of Algorithm 3.

Numerical results

We now present results obtained with a code (3 dimensionnal in velocity) developped following methods proposed in Section 2 to solve model [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF] with initial condition (19). The initial density functions f ini D (v) and f ini T (v) are Maxwellian distributions of density ρ ini a , velocity u ini a and temperature θ ini a that is to say

f ini a (v) = ρ ini a (2πθ ini a ) 3/2 exp - |v -u ini a | 2 2θ ini a , a ∈ {D, T }. (46)
We use the physical values given in Table 1 and the values of numerical parameters given in Table 3. Moreover, the initial density functions f ini α (v) and f ini n (v) are equal to zero. The initial time step ∆t 0 is chosen in such a way condition (33) is satisfied for the proposed test-cases.

m D m T m α m n Q (kg) (kg) (kg) (kg) (MeV) 3
ρ ini D ρ ini T θ ini D θ ini T u ini D u ini T (keV) (keV) (m•s -1 ) (m•s -1 )
10 20 10 20 15 20 0 0

N 0 D ∆t 0 (s) 10 4 10 -1
Table 2: Initial values of macroscopic quantities and initial numerical parameters.

Validation with a Maxwellian cross-section

We first validate the numerical method proposed in §2.2 for a Maxwellian cross-section [START_REF] Lapeyre | Méthodes de Monte-Carlo pour les équations de transport et de diffusion[END_REF], with λ = 10 -21 m 3 •s -1 , thanks to the explicit time evolution of macroscopic quantities obtained in §1. [START_REF] Duane | Fusion cross-section theory[END_REF]. We can see on Figure 2 that the evolution of the number densities of all species are very close to the one of the explicit solution (24).

f h,n a,1 (t n , v) = (S 3,h * f n a,1 )(v) = ˆR3 S 3,h (v -v ′ )f n a,1 (v ′ )dv ′ = N n a i=1 ωS 3 (v -V n i,a ),
where

S 3,h (x) = 1 6h            2 - x h 3 if h ≤ |x| < 2h, 4 -6 x h 2 + 3 x h 3 if 0 ≤ |x| < h, 0 otherwise.
We can observe that if initialy f a,x (a ∈ {D, T }) is a Maxwellian distribution, it seems that it remains a Maxwellian distributions at t > 0 which is in agrement with (20). Moreover, Figure [START_REF] Duane | Fusion cross-section theory[END_REF] shows that

Simulation with the fusion cross-section

We consider here the model ( 1) where the cross-section is a fusion cross-section σ f given by ( 8) with the values of Duane coefficients given by Table 3. Figure 5 presents the time evolution of the Figures 7 and8 present the reconstructed density functions respectively of the first component of the velocity and of the energy of each species. Figure 7 shows that the density functions in velocity of α and n are clearly not Maxwellian distributions. Figure 8 shows that the energy density functions of α and n seem quite close to a Maxwellian distribution and are centered around the standard values (16).

Finally, we compare on Figure 9 the computed fusion reaction rate τ with explicit expression (23) which corresponds to a Maxwellian cross-section. Since the parameter λ in (23) is not defined for fusion cross-section (23), we compute it from the value τ 0 fusion of τ computed at t = 0 for the fusion cross-section: thus, we set λ = τ 0 fusion /(ρ ini ) 2 . This explicit expression gives a quite good approximation of the computed values of τ as we see on Figure 9. This expression should be compared to results of computations of the fusion rate obtained in [START_REF] Goncharov | Practical Calculation of Nuclear Fusion Power for a Toroidal Plasma Device with Magnetic Confinement[END_REF] or in [START_REF] Heltemes | Analysis of an Improved Fusion Reaction Rate Model for Use in Fusion Plasma Simulations[END_REF].

On the opportunity to use the asymptotic model

We now discuss the opportunity to use the asymptotic model (18) instead of model [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF]. Figure 10 presents the reconstructed energy density function of α particles for different values of θ = min(θ D , θ T ) at time t = 10 7 s (at this time, almost every possible fusion reactions have been made and f α (t, v) is close to lim t→+∞ f α (t, v)). Initial conditions are again given by (19) with the initial density f ini D and f ini T given by (46) with parameters (3) except for the initial kinetic temperatures θ ini D and θ ini T . Indeed, we now take θ ini D = θ ini T ∈ {5, 10, 15, 20} keV. We see that the sprawl of the distributions in energy of α particles becomes smaller when θ decreases. This confirms at the discrete level that the energy density functions of α becomes closer to a Dirac mass centered at energy E α (given by ( 16)) when the parameter ε defined by (15) tends to zero. We verify that this numerical result is also valid for the neutrons (by replacing E α with E n ). Consequently, when the temperature of the plasma (that is to say θ ini D and θ ini T ) is not too high, it seems to be valid to approach model (1) by (18). However, as explained in Remark 6, the numerical cost of the resolution of model ( 18) is the same as the one of the resolution of model (1) because of the calculus of the fusion reaction rate τ (t). Then, an interesting possibility to compute evolution of α particles and neutrons could be to use (45) with the explicit expression of (23) of τ (t) since Figure 9 shows that this expression -with an adequate value of λ -is a good approximation of the fusion reaction rate. Nevertheless, we should make further investigations to compare the evolution of the distributions in velocity of α particles and neutrons of model (45) with model (1), in particular in the case of non Maxwellian initial Deuterium and Tritium distribution densities.

Conclusion

We have proposed and justified a Monte-Carlo algorithm -adapted from a Monte-Carlo algorithm for the Boltzmann equation -to resolve the homogeneous kinetic model (1) describing the nuclear fusion reaction between a Deuterium ion and a Tritium ion and giving an α particle and a neutron. We have also developed a spacially homogeneous code (3D in velocity) based on this Monte-Carlo algorithm. We can thus obtain the evolution of the distributions in velocity, distributions in energy, and macroscopic quantities of Deuterium and Tritium ions, α particles and neutrons starting from initial density functions. The initial distributions used in the presented examples are only Maxwellian distributions but we can also use other initial distributions such as a sum of Maxwellian distributions with different densities, velocities and temperatures.

We validated the proposed Monte-Carlo algorithm thanks to the use of a Maxwellian cross-section instead of the fusion cross-section. Indeed, the use of a Maxwellian cross-section allows to obtain explicit solutions. Moreover, the numerical cost of our computations is very moderate: a few minutes for a computation with physical and numerical parameters of Table 1 and pratically immediate if the number of initial numerical particles N 0 D is 10 3 instead of 10 4 . However, it could be interesting to have a model less costly in order to use it in non-homogeneous simulations with thousands of cells. If we are only interested in the evolution of α particles -because they are the ones who can trigger instabilities in the plasma -and neutrons, a possibility could be to use asymptotic model (45) with an explicit value of the fusion reaction rate τ (t). This method is already used as for example in [START_REF] Farengo | Simulation of alpha particle current drive and heating in spherical tokamaks[END_REF]. Here, we proposed to approach the fusion reaction rate τ (t) by using the explicit formula obtained when the cross-section is a Maxwellian cross-section. Nevertheless, further investigations should be made to validate this approach, especially when the plasma is not initially at the thermodynamical equilibrium.

Figure 2 :

 2 Figure 2: Time evolution of number densities ρ a , a ∈ {D, T, α, n}: comparison between the numerical solutions obtained by Algorithm 2 and explicit solutions given by expressions (24).

Figure 3

 3 Figure 3 presents the density functions f h D,1 and f h T,1 of the first component of the velocity of D and T ions at different times. The density functions are reconstruted from the values of the velocities of numerical particles with a convolution with a B-splines of order 3 like in [15]:

Figure 3 :

 3 Figure 3: Reconstructed density functions of the first component of the velocity of Deuterium and Tritium ions at different instants in the case of the Maxwellian cross-section.the kinetic temperatures of every species are constant as it is predicted by (26).

Figure 4 :

 4 Figure 4: Time evolution of kinetic temperatures in the case of the Maxwellian cross-section.

Table 3 :

 3 • 10 -24 1.368 • 10 -2 1.076 4.09 • 10 -26 Values of Duane coefficients in fusion cross-section σ f (8). number densities of D and T ions on one hand and of α particles and neutrons on the other hand, and the evolution of the total energy of the plasma E tot .

Figure 5 :

 5 Figure 5: Time evolution of number densities ρ a , a ∈ {D, T, α, n} and of the total energy in the case of the fusion cross-section.

Figure 6

 6 Figure6presents the time evolution of kinetic temperatures (in keV) of the species. We observe that the kinetic temperatures of D and T are not constant but are decreasing, contrary to the case of the Maxwellian cross-section (compare with Figures 4). Moreover, if we compute the values of the kinetic energy of α particles and of neutrons from kinetic temperatures T α and T n of Figure6, we

  Figure6presents the time evolution of kinetic temperatures (in keV) of the species. We observe that the kinetic temperatures of D and T are not constant but are decreasing, contrary to the case of the Maxwellian cross-section (compare with Figures 4). Moreover, if we compute the values of the kinetic energy of α particles and of neutrons from kinetic temperatures T α and T n of Figure6, we

Figure 6 :

 6 Figure 6: Time evolution of kinetic temperatures in the case of the fusion cross-section.

Figure 7 :

 7 Figure 7: Reconstruted density functions of the first component of the velocity for each species in the case of the fusion cross-section.

Figure 8 :

 8 Figure 8: Reconstruted energy density functions in the case of the fusion cross-section.

Figure 9 :

 9 Figure 9: Comparison between the time evolution of the computed fusion rate τ in the case of fusion cross-section and the explicit expression (23) of τ in the case of the Maxwellian cross-section.

Finally, we

  wish to underline that this study is a first step for the study of a more complet model with feedback effect of α particles on Deuterium-Tritium plasma.

Figure 10 := λ max ∆tω 2

 102 Figure 10: Reconstructed energy density function of α particles for different initial values of θ = min(θ D , θ T ).

Table 1 :

 1 .34 • 10 -27 5.01 • 10 -27 6.68 • 10 -27 1.67 • 10 -27 17.59 Physical constants.

fa(t, va) is in s 3 •m -3

It is also possible to obtain an expression of Gα on a strong form[START_REF] Dellacherie | Nuclear Collisions Models with Boltzmann Operators[END_REF] but we do not need it here.

Other expressions of fusion cross-section could be found, such the one of[START_REF] Li | Fusion cross-sections for inertial fusion energy[END_REF].

The Maxwellian cross-section, which corresponds to the cross-section of the model of "Maxwellian molecules" for Boltzman operator, should be considered as a theoretical model. However, it is widely used in kinetic theory (see[START_REF] Villani | A review of mathematical topics in collisional kinetic theory[END_REF]) to perform many explicit calculations.

Acknowledgment

The autors would like to thank the financial support of the Large Scale Initiative FUSION of INRIA.

Annex A: proof of Propositions 3 and 4

Proof. We use the fact that for every integer N ≤ min(N n D , N n T ) and every function g : {1, . . . , N n D } × {1, . . . , N n T } → R, we have

where P is the set of possible ways to make N distinct pairs (r k , s k ) {1≤k≤N } among {1, . . . , N n D } × {1, . . . , N n T }. For each test function ϕ, we get ∆t

where N n is defined by (35) and where we have set


J Delahaie

1

Sociopragmatic competence in FFL language teaching : discourse markers as contextualisation cues

The common European Framework of Reference for Languages (1999) : focus on pragmatics and sociolinguistics skills, i.e. :

-knowledge of scripts, micro-and macro-functions -ability to manage an interaction -knowledge and use of relational and social markers Neither pragmatic nor sociolinguistic skills are independent of the linguistic components of communicative competences : certain linguistic features convey information about context (the extralinguistic environment of the utterance). The aim of this presentation is to demonstrate that discourse markers are contextualisation cues and to propose some ways in which discourse markers (DM) might therefore be taught. Data : authentic oral interactions from the Lancom corpus (Université of Louvain-Leuven), available on line at : http://bach.arts.kuleuven.be/elicop/ Transcriptions of 13 interactions in a travel agency in Marcq-en-Baroeul (Lille), 30.000 words (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006).

I. Some definitions : sociopragmatic competence, discourse markers I.1. Sociopragmatic Competence [START_REF] Leech | Principles of pragmatics[END_REF], differenciates :

-sociopragmatic knowledge : the contextual and social variables that determine appropriate pragmalinguistic choices.

-pragmalinguistic knowledge : particular linguistic resources for formulating a particular type of speech act. Rose (199 : 171) : "The ability to use available linguistic resources (pragmalinguistics) in a contextually appropriate fashion (sociopragmatics)."

I.2. Discourse markers as contextualisation cues

Discourse markers are : -words that have a pragmatic function and which are not grammatically integrated into a sentence -linguistic expressions that act as various parts of speech (adverbs, conjunctions, interjections) -words that perform principally expressive and social functions : they convey information about speakers' attitudes and actions, about shared knowledge (or lack of it), about the exchange structure, and so ons.

Various theories deal with discourse markers :

-"discourse analysis", "data-driven approach" of [START_REF] Schiffrin | Discourse markers[END_REF][START_REF] Schiffrin | Discourse Markers : Language, Meaning, and Context[END_REF] : data is used as a starting point for an account of the use and distribution of markers in everyday discourse, and to outline a pattern of use.

-"theory driven-approach" : pragmatics integrated into semantics and therefore the fundamental meaning of lexical items (especially Anscombre) (a non truth-conditional theory of language which suggests that to speak is first and foremost to have an effect upon another person).

II. How should sociopragmatic competence and discourse markers be taught (and learned)?

Sociopragmatic competence [START_REF] Bardovi-Harlig | Empirical evidence of the need for instruction in pragmatics[END_REF] demonstrates the benefits of instruction versus immersion in various aspects of pragmatics. But how might pragmatic skills be taught? -enhance learners "metapragmatic consciousness", invite the teacher to explain what remains implicit in the target language [START_REF] Ball | Colloquial French grammar. A practical guide[END_REF][START_REF] Bardovi-Harlig | Empirical evidence of the need for instruction in pragmatics[END_REF] -make discourse markers relevant : Kasper and Roever (2005 : 318), in order to be acquisitionally available, the input must seem relevant and worthwhile. Learners' attention should be drawn "to the linguistic forms by which [an] action is implemented, its immediate interactional or textual context."

Discourse markers

A problem : the multifunctionality of DM DM are appropriate in various contexts whose parameters are too numerous to be demonstrated practically. Illustration : the DM voilà. This very frequent DM in French has numerous functions, and might index various contextual parameters : -as an affirmative response to a request for confirmation :

(1). E-c'était l'hôtel X hein c'est ça hein ? C-voilà à Lesbos et il y a une semaine à Kyos (travel agency, « Lesbos », 2006) Such confirmation requests should not be confused with genuine questions, to which voilà could not serve as an answer :

(2). -est-ce que c'est l'hôtel X ? -oui -*voilà/*d'accord -its second most common function is to conclude, at the end of a turn in dialogue :

(3) E-[…] alors on va regarder ce qu'on peut ce que je peux regarder parce que en fait si vous voulez le samedi pendant les vacances les euh touropérateurs ferment plus tôt (nom) ferme à seize heures trente et (nom) ferme à dix-sept heures et euh voilà (travel agency, « Chypre », 2006) This occurs especially frequently during what we might call "sequences of clusure" of an interaction :

(4). C-bon j/e peux m/e permettre de l'emmener ? E-bien sûr sans souci sans souci C-à c/e moment-là j/e vais le j/e vais cogiter un peu E-et donc euh voilà sachez qu/e pour l'aérien donc i/l y a pas de souci C-oui oui E-en aller simple en aller-retour ou en retour simple C-tout à fait E-et qu'à côté d/e ça euh si c'est pour l'hiver effectivement euh ce s/era plus eu:h à partir du mois d/e septembre C-pour l'hiver oui d'accord E-OK ? C-donc j/e verrai demain je n/e sais pas j/e vais voilà (travel agency, « Suède », 2006)

-more rarely, voilà is used as an overture at the beginning of the exchange : a means of opening a conversation :

(5). [beginning of a telephone call] E-alors en fait voilà je voulais vérifier parce que j'aurais éventuellement une personne qui souhaiterait se rajouter […] (travel agency, « Lesbos », 2006)

How should these usages be taught? -make lists of DM?

The danger is the "paradigmatic curse" identified by [START_REF] Beacco | Un rendez-vous manqué? Théorie du discours et grammaire en didactique[END_REF] concerning vocabulary in foreign language teaching : giving learners lists of words to memorize, completely out of context. This is nevertheless precisely what the Référentiels pour le français do : provide lists of sentences or DM, linked to particular speech actsfor example to "confirm" : oui, tout à fait, si, bien sûr que si ; or to "conclude" : (bon ben) voilà quoi, voilà, enfin, donc, bref, etc. (Niveau B2 pour le français, un référentiel, 2004 : 72 and 114).

-Give a semantic or syntactic explanation ? The problem is that DM like voilà have neither a specific denotation nor a clear syntactic role. They have, instead, a subjective/expressive meaning [START_REF] Schiffrin | Discourse Markers : Language, Meaning, and Context[END_REF] : they convey information about the speakers' relationship, their shared knowledge (1 st use of voilà), and take part in the co-construction of an interaction (2 nd use of voilà).

-Increase awareness? Some researchers suggest increasing learners' awareness of DM [START_REF] Wichmann | Discourse markers : a challenge for linguists and teachers[END_REF], [START_REF] Dewaele | L'acquisition de la compétence sociopragmatique en langue étrangère[END_REF] : teaching with authentic spoken data (video, transcriptions). This all seem to support the idea that DM cannot be explained and are therefore not teachable. III. Teaching DM in context : other propositions (the contrastive approach, and models of exchange and dialogue) III.1. Theoretical starting point DM have very complex meanings, but not (any longer) denotations : they convey information about the context. Schiffrin's data-driven theory (2001 ;[START_REF] Schiffrin | Discourse markers[END_REF] : DM are contextualisation cues, which act upon different aspects of an interaction (cognitive, expressive, social, and textual levels). Markers are therefore multifunctional. How can teachers explain the effects of multiple contextual constraints on DM?

Another theoretical starting point : DM maintain a core of semantic traits throughout (and despite) their multiple functions. They convey consistent contextual indices of speaker relationships and the construction of an interaction. This is the theoretical point of view that will be adopted here to examine the multifunctional French DM voilà. In particular, the presentation will employ Anscombre's framework of Integrated Pragmatics and Stereotype Theory.

III.2. A contrastive approach to the discourse marker : the example of voilà

Our hypothesis is that functional variations of voilà are surface effects of the context, but voilà has a unique semantic meaning : Voilà confirms or concludes the interlocutor's utterance, marking it as conforming (as it were) to a stereotype ; it is something that is expected and "follows" (thereby representing a kind of shared knowledge) ; it is a very strong consensual marker.

N.B : Stereotype theory : stereotypes are an "open series of utterances attached to a lexical item, which constitute its meaning" (Anscombre, 1995 : 82) ; according to this theory, speaking necessarily involves the implication of one or more stereotypical utterances, which are culturally determined, and which presuppose certain types of sequences and conclusions. Stereotype theory exists within the framework of Integrated Pragmatic Theory whereby the semantic meaning associated with a sentence is said to include indications of a "set of instructions related to the way in which the utterance is to take its place within a discourse dynamic" (Anscombre, 1990 : 50). Despite the fact that it is not normally used to explain interactive aspect of linguistics, the interest of this theory for the present study is to be found in its focus upon sequences of utterances, and in its hypothesis that these sequences are often based upon the presuppositions carried by words in the form of "stereotypes".

Core meaning of voilà

When voilà confirms : consensual meaning based on speakers' shared knowledge. When voilà concludes : it presents the situation as part of a stereotypical and expected sequence, sometimes in order to replace an explicit explanation. When voilà is employed within a monological conclusion, a sentence can be finished without proceeding to a normal grammatical closure :

(6) E-[…] alors on va regarder c :e qu'on peut c/e que j/e peux regarder parce qu/en fait si vous voulez le samedi pendant les vacances les eu :h tour-opérateurs ferment plus tôt X ferme à seize heures trente et Y ferme à dix-sept heures et euh voilà C-d'accord mais déjà euh qu'est-c/e qu'on pourrait en fait on est jamais parti avec des tour op […] (travel agency, "Chypre", 2006)

In the travel agency interactions, when voilà appears at the end of the interaction, it is used at the end of a script (with a set of expected sequences, exchanges and speech acts) as a conclusion : a usage that is less frequent in conversational interactions.

Application 1 : contrastive exchanges

In ordre to highlight the relevance of the semantic and pragmatic meaning of voilà, contrastive, binary models of exchanges might be proposed that focus on two different discourse markers. For example : voilà and d'accord as a binary pair of responsive speech acts : Application 2 : models of contrastive interactions How might voilà be presented as a concluding or opening marker? A contrastive approach through dialogue : voilà, as a consensual marker, is used at the end of a script , as a confirmation of an expected conclusive sequence. It rarely occurs in conversational interactions (without a script), replaced by tiens (opening) or bon (conclusion), as Traverso has demonstrated (1996). Thus a teacher might want to compare two different types of interaction : conversation familière : rencontre en rue de copains (interaction de Lancom remaniée) interaction à script : à l'agence de voyage (interaction de Lancom remaniée) (C1 et C2 sont dans la rue et se rencontrent) C1-tiens salut Bertrand, comment ça va ? C2-très bien et toi ? C1-bah ça va je suis en train de faire de courses là C2-tiens au fait tu as des nouvelles de Jérémy il devait partir en Angleterre C1-oui il a réussi à partir pour quinze jours, il rentre dans une semaine C2-ah oui c'est bien C1-bon ben il faut que j'y aille de toute façon on s'appelle hein ? C2-d'accord allez salut !

Conclusion

Teaching DM in FFL is impossible without a systematic linguistic description of the DM. This description should demonstrate that DM prescribe linguistic systems and do not vary according to external contexts : they explain contexts. The project of the Dictionnaire raisonné des marqueurs discursifs du français, undertaken by D. Paillard (see [START_REF] Paillard | Prise en charge, commitment ou scène énonciative[END_REF], in collaboration with the University of Hanoï and the Agence universitaire de la Francophonie : to describe and explain DM in French as a linguistic and pragmatic system (according to Culioli's Enunciation Theory), is to be followed by a thorough didactic model for teaching and learning of this system.

  (7) . C-là vous attendez une confirmation d'AQUATOUR ? E-voilà c-d'accord (travel agency, « un père organisé », 1993) D'accord is frequently the third "evaluative" move of an exchange (36,9% of the evaluative responses in our data), it serves as a ratification of new information.Dans une agence de voyage, C est le client, et E l-ce que vous avez une idée de voyage ? C-oui Figure1. Des modèles d'échange contrastifs

(

  Figure 2. Des modèles d'interaction contrastifs

i-mean@uwe2 Conference, atelier AFLS, 13-15 April 2011, Bristol