S Brull

P Charrier
email: pierre.charrier@math.u-bordeaux.fr

L Mieussens
email: luc.mieussens@math.u-bordeaux.fr

Nanoscale roughness effect on Maxwell-like boundary conditions Nanoscale roughness effect on Maxwell-like boundary conditions for the Boltzmann equation

Keywords: numbers: 51.10, 34.35, 02.70 gas-surface interaction, kinetic theory

It is well known that the roughness of the wall has an effect on microscale gas flows. This effect can be shown for large Knudsen numbers by using a numerical solution of the Boltzmann equation. However, when the wall is rough at a nanometric scale, it is necessary to use a very small mesh size which is much too expansive. An alternative approach is to incorporate the roughness effect in the scattering kernel of the boundary condition, such as the Maxwell-like kernel introduced by the authors in a previous paper. Here, we explain how this boundary condition can be implemented in a Discrete Velocity approximation of the Boltzmann equation. Moreover, the influence of the roughness is shown by computing the structure scattering pattern of mono-energetic beams of the incident gas molecules. The effect of the angle of incidence of these molecules, of their mass, and of the morphology of the wall is investigated and discussed in a simplified two-dimensional configuration. The effect of the azimuthal angle of the incident beams is shown for a three-dimensional configuration. Finally, the case of non-elastic scattering is considered. All these results suggest that our approach is a promising way to incorporate enough physics of gas-surface interaction, at a reasonable computing cost, to improve kinetic simulations of micro and nano-flows.

I. INTRODUCTION

In some applications we have to consider flows at high Knudsen number on domains bounded by walls. Examples of such situations are external flows at high altitude around re-entry bodies or internal flows in micro or nano devices (MEMS or NEMS) and in sub-micrometer porous media. In such configurations the flow is well described by a kinetic equation (the Boltzmann equation or the BGK model), together with boundary conditions imposed on the wall. With more and more powerful computers such models are now currently used for numerical simulations. A boundary condition for the Boltzmann equation is a way to take into account, at the kinetic scale, the complex gas-surface interaction which takes place at a very small scale (typically on less than one nanometer). The most used boundary conditions are the Maxwell and the Cercignani-Lampis conditions (see [START_REF] Cercignani | The Boltzmann Equation and its Applications[END_REF] for a detailed review on boundary conditions and 2 3 for recent works). However none of them take into account some small or medium scale roughness of the wall, though it is well-known that such roughness may notably affect the flow near the wall. As a result, those standard boundary conditions may turn out to be unable to describe with enough accuracy the transport properties of the flow near the surface.

To overcome this difficulty for the Maxwell condition, a mathematical framework, recently introduced by the authors [START_REF] Brull | Gas-surface interaction and boundary conditions for the Boltzman equation[END_REF] , allows to derive boundary conditions for a periodic rough wall by using a rigorous upscaling method. In the simplest case, corresponding to a flat wall, this condition looks like a Maxwell condition but with accommodation coefficients that depend on the velocity of the impinging molecules and on the temperature and morphology of the wall (see another recent approach 2 that also leads to variable accommodation coefficients). At the small scale the gassurface interactions are assumed to be governed by a potential energy surface (PES) which was not specified in [START_REF] Brull | Gas-surface interaction and boundary conditions for the Boltzman equation[END_REF] .

In the present paper we propose a simple but physically reasonable PES and a way to incorporate it in this mathematical framework, in order to get an improved boundary condition for the Boltzmann equation easy to implement for numerical simulations with a discrete velocity method (DVM). Many contributions on the gas-surface interaction have been published in past years either using very simple models, like the "hard cube" model [START_REF] Logan | Simple Classical Model for the Scattering of Gas Atoms from a Solid Surface[END_REF] or the washboard model [START_REF] Tully | Washboard model of gas-surface scattering[END_REF] or much more elaborate ones 7 8 . Here we suggest to use an intermediate model, simple enough to be tractable, but precise enough to include at the kinetic scale the effect of van der Waals forces and of roughness at atomic scales. The PES is defined by a lattice of Lennard-Jones atoms. In an auxiliary problem, for given incident velocities of impinging molecules, we compute the trajectories of a nanoscopic sample of gas atoms in a periodic cell of the surface layer, in which gas molecules interact with surface atoms by the simplified PES. From those data computed once and for all, we get an approximation of the reflection kernel of the new boundary condition suitable for DVM numerical simulations of the Boltzmann equation. At the same time we obtain the scattering pattern of mono-energetic beams of incident gas atoms by the wall and we can visualize and analyze the ability of our approach to incorporate at the kinetic level information from the atomic scale.

The paper is organized as follows: Section II reminds the modified Maxwell-like boundary condition 4 (MLBC) and its derivation for elastic scattering and non-elastic scattering as well. In section III we describe how to implement this MLBC in a DVM context. Section IV is devoted to the numerical computation of scattering patterns to analyze the effect of nanoscale roughness. The impact of the key parameters (incident angle and energy of the impinging gas atoms, parameters of the Lennard-Jones potentials) on the so-called structure scattering pattern (elastic scattering) is shown (for 2D and 3D configurations) and some indications are given for non-elastic scattering. The paper closes with some comments and conclusions in section V.

II. NEW BOUNDARY CONDITIONS FOR THE BOLTZMANN EQUATION

A. A nanoscale gas-surface interaction model

We start from a nanoscale model describing the gas-surface interaction which generalizes to rough walls the model introduced in [START_REF] Borman | Theory of nonequilibrium phenomena at a gas-solid interface[END_REF] and used in [START_REF] Aoki | A hierarchy of models related to nanoflows and surface diffusion[END_REF] . In this model the interaction between the wall and the single atom gas molecules through van der Waals forces are taken into account in a thin surface layer (with thickness L typically smaller than one nanometer). In all the following (except at the end of section IV A), for the sake of simplicity, we assume that the atoms move in a 2D half-plane and we consider the following configuration: the solid occupies the half-space z > L, the gas phase is constituted by the gas atoms in the half-space z < 0, outside of the range of the surface forces, and we consider separately the surface layer 0 < z < L, where the gas atoms move within the range of the surface potential. The gas flow in this surface layer is modeled by the Boltzmann equation (or a BGK-type model equation), with a Vlasov term to take into account the part of the interaction that depends on the frozen position of the atoms of the solid wall (the long range interactions), and a collision term between gas atoms and phonons to take into account the thermal fluctuations of the atoms of the solid (short range interactions) (see [START_REF] Borman | Theory of nonequilibrium phenomena at a gas-solid interface[END_REF] for a physical justification of this approach). We consider (as in 4 ) the following general configuration for the wall: the interaction potential V(x, z) is such that

V(x, z) = V # ( x L * , z), (1)
where L * = L λ * and λ * is a positive constant that characterizes the roughness of the wall, and V # (y, z) is a periodic function of the nanoscopic variable y with period 1. This nanoscopic variable y allows us to describe how a gas atom impinging the surface layer at a microscopic coordinate x sees the nanoscopic roughness of the wall. Moreover we assume that there exists a 1-periodic function

z = ζ ∞ (y) with 0 ≤ ζ ∞ (y) < L such that lim z<ζ∞(y),z→ζ∞(y) V # (y, z) = +∞.
(2)

V # (y, z) = V m , for z ≤ 0. (3)
Condition ( 2) is a mathematical way to ensure that no gas atom can penetrate the wall, i.e. that it is a non porous wall and condition (3) means that the potential is constant outside the surface layer (see Figure 1). With these assumptions, the flow

z 0 < y < 1 V # (y, ζ ∞ (y)) = +∞ Surface layer z = L z = 0 V # (y, z) = V m z = ζ ∞ (y) Solid FIG.
1. Parametrization of the surface layer of gas atoms is described by the following system of kinetic equations

∂ t f + v.∇ x,z f = Q m,m [f, f ], z < 0, ∂ t f + v.∇ x,z f - 1 m ∇ x,z V(x, z).∇ v f = Q m,m [f, f ] + Q m,p [f ], 0 < z, (4)
where Q m,m is the Boltzmann collision operator (or a BGK like relaxation operator)

and where the gas atom-phonon collision term reads (see for instance [START_REF] Yu Krylov | One dimensionnal surface diffusion. II Density dependence in a corrugated potential[END_REF] )

Q m,p [f ] = K(E(x, z, v), E(x, z, v ); T w )(M w (E(x, z, v))f (v ) -M w (E(x, z, v))f (v))dv , (5)
with

M w (E(x, z, v)) = exp - m|v| 2 2kT w - V(x, z) kT w , (6)
and K(E, E ; T w ) is the probability of transition per unit time from the state E to the state E in a "collision" with a phonon. Moreover Q m,p satisfies the usual property of conservation of mass, an H-theorem and a symmetry property, and it has a Maxwellian equilibrium. Since we assume in all the following that the wall temperature T w is known and fixed, we will drop the dependence of K on T w .

Finally we introduce two times, τ ms and τ f l , that will play a key role in the following and that are defined as follows: τ ms is the time of relaxation of gas atoms by the phonons and is defined by

τ ms = τ ms (E(x, z, v)) = K(E, E )M w (E )dv -1 , (7)
where E = E(x, z, v ), while τ f l (v) is the time of flight of a gas atom across the surface layer (without collision with a phonon): it depends on the incident velocity of the gas atom at the entrance of the surface layer, on the potential, and on the morphology of the wall.

B. Derivation of boundary conditions for elastic scattering on a rough wall

The derivation of boundary conditions for the Boltzmann equation on a wall with nanoscale roughness relies on a two scale asymptotic analysis from the above gassurface interaction model. In this section we only indicate the main ideas of this approach, skipping the technical mathematical details. The reader interested in a more rigorous presentation can refer to ( 4 ).

In many applications we have to consider rarefied flows on microtubes with a typical length of 1 to 10 cm (see for instance [START_REF] Perrier | Mass flow rate measurements in microtubes: from hydrodynamic to near free molecular regime[END_REF] ). In such configurations a reasonable mesh size l * for a simulation based on the Boltzmann equation would be around 10 µm, so that L * l * . Associated with those two space scales are two time scales t * B , the reference time of evolution of the Boltzmann equation defined by t * B = l * /v * and τ * f l , the reference time of flight of a gas atom across the surface layer defined by τ * f l = L * /v * where v * is a typical velocity of atoms of the gas flow (for instance v * = kT w /m). Then we have τ * f l /t * B = L * /l * 1. In most applications we can also consider that the reference relaxation time τ * ms (that can be defined as τ ms (E(x, 0, v * ))), is such that τ * ms t * B , but the interaction between the wall and the gas atoms depends mainly on the respective order of magnitude of the characteristic time scales τ * f l and τ * ms . In this subsection we assume that τ * f l τ * ms , that is to say that the time of flight of a gas atom across the surface layer is small enough so that the probability to have a collision with a phonon can be neglected. This is a limit regime but it can be relevant for low wall temperature or for light gas molecules.

Since the potential V is rapidly oscillating in the x-direction at the nanoscale, so is the solution of the system (4), at least near the surface layer. Thus a usual technique in the homogenization method is to look for a solution of (4) in the form (see 13 14 )

f (t, x, z, v) = F (t, x, x/L * , z, v), (8)
where the function F (t, x, y, z, v) is periodic with respect to y with period 1. Note that with such a definition, we have

∂ x f = ∂ x F + 1 L * ∂ y F, ∂ x V = ∂ x V # + 1 L * ∂ y V # . (9)
In the surface layer (for z > 0), the variations of the function F with respect to the variables y and z are significant while we can neglect the (much slower) variations of F with respect to t and x. And the right-hand side can also be neglected because we neglected the probability of collision between gas atoms and phonons (and even more the probability of collision between gas atoms). In the bulk flow we are interested in describing the variation of F at the microscale (i.e. l * ) and the fast variation at the nanoscale (L * ) can be taken into account only in average so that the bulk flow is described by the following distribution function

g(t, x, z, v) = 1 0 F (t, x, y, z, v) dy, for z < 0. ( 10
)
Thus, in the limit regime we consider in this section, the gas flow can be described by the following simplified sytem

∂ t g + v x ∂ x g + v z ∂ z g = Q m,m [g, g], z < 0, (11)
g(t, x, 0, v x , v z ) vz<0 = 1 0 F (t, x, y, 0, v x , v z )dy, (12)
v x ∂ y F + v z ∂ z F - 1 2 ∂ y V # (y, z)∂ vx F - 1 2 ∂ z V # (y, z)∂ vz F = 0, 0 < z, (13)
F (t, x, y, 0, v x , v z ) |vz>0 = g(t, x, 0, v x , v z ), ∀y ∈ [0, 1]. ( 14
)
This coupled system must be understood as follows. Equation ( 13) defines the distribution F of gas atoms in a periodic cell of the surface layer: it depends on the distribution g(t, x, 0, v x , v z )| vz>0 of gas atoms coming from the bulk flow that enter into the surface layer with v z > 0. Then the distribution F (t, x, 0, v x , v z )| vz<0 of gas atoms that leave the surface layer (with v z < 0) can be, in principle, computed by solving ( 13)-( 14). Finally, the distribution g of gas atoms in the bulk flow is defined through ( 11)- (12). Relation (12) means that the distribution of gas atoms that come into the bulk flow from the surface layer is the average on one periodic cell of the distribution F (t, x, 0, v x , v z )| vz<0 of gas atoms that leave the surface layer. In other words, relation (12) can be viewed as an implicit boundary condition for the bulk flow governed by (11), which requires to solve the boundary layer problem ( 13)-( 14). This is summarized in Figure 2. Now, we explain how a boundary condition for the bulk flow can be constructed. In fact, (13) can be written d ds F (t, x, y(s), z(s), v(s)) = 0 ( where s → (y(s), z(s), v(s)) is the trajectory of a gas atom across the surface layer, see equations ( 24)-(25) in section III). Then ( 13)-( 14) can be integrated exactly, and the outgoing value F (t, x, y, 0, v x , v z )| vz<0 can be written as a function of the incoming value (14). This outgoing value can be used in the implicit boundary condition (12), that can now be written explicitly as

g(t, x, 0, v) |vz<0 = 1 0 g(t, x, 0, -Λ 2 (y, -v)) dy, (15)
where Λ(y

, v ) = (y = Λ 1 (y , v ), v = Λ 2 (y , v ))
is the mapping that gives the position and velocity (y, v) of a gas atom that leaves the surface layer (i.e. with

z z = 0 z = L eq. (11) for g(t, x, z, v x , v z ) g(t, x, 0, v x , v z )| vz>0 g(t, x, 0, v x , v z )| vz<0
eq. ( 13) for F (t, x, y, z, v x , v z )

F (t, x, y, 0, v x , v z )| vz>0 periodic cell
Surface layer Solid eq. ( 14)

F (t, x, y, 0, v x , v z )| vz<0
eq. ( 12)

FIG. 2. Links between the bulk flow equation ( 11) and the boundary layer problem ( 13)-

v z < 0 at z = 0) as a function of its position and velocity (y , v ) when it enters into the surface layer (i.e. with v z > 0 at z = 0), see figure 3. This mapping is defined through the equations of the trajectory of a gas atom (see section III and equations ( 24)-( 25)). In other words, with a probabilistic interpretation, the incoming condition (14) means that a gas atom impinging the surface layer at point (x, z = 0) (at microscale), sees the nanoscale roughness from point (x, y, z = 0), y ∈ [0, 1] with uniform probability in y. Equation (13) gives the deterministic trajectories of the gas atoms inside the surface layer under the action of the interaction potential. Finally the outgoing distribution function at microscale is defined as the mean (in y) of the outgoing nanoscale distribution functions at point (x, y, z = 0) (14). After some algebra the reflection kernel of condition ( 15) is found to be

z = L z = 0 (y ′ , v ′ ) = (Λ 1 (y, -v), -Λ 2 (y, -v)) with v ′ z > 0 with v ′ z > 0 (y ′ , -v ′ ) (y, -v) (y, v) = (Λ 1 (y ′ , v ′ ), Λ 2 (y ′ , v ′ )) y = 1 y = 0
g(t, x, 0, v)|v z | |vz<0 = v z >0 R(v → v)g(t, x, 0, v )|v z | dv , (16)
where

R(v → v) = |v z | |v z | 1 0 δ(v + Λ 2 (y, -v))dy, (17)
and δ(y) is the Dirac delta function. This condition satisfies the properties of nonnegativeness, normalization and reciprocity and also a H-theorem, with an entropy production vanishing only for a flat wall (see [START_REF] Cercignani | The Boltzmann Equation and its Applications[END_REF] for a definition of these notions).

Formula ( 17) clearly shows that the scattering kernel R(v → v) essentially depends on the trajectories of the gas atoms inside the surface layer (through the mapping Λ).

See section III for an application of this reflection kernel with a given interaction potential. Note that it can be easily checked that, for a flat wall (i.e. for an interaction potential V(x, z) = V(z)), the boundary conditions ( 16) reduces to the classical specular reflection.

C. Derivation of boundary conditions for "non-elastic" scattering on a rough wall

We consider non-elastic scattering and in a first step we assume that

τ * ms τ * f l t * B . (18)
This means that during their flight in the surface layer, some gas atoms have time to have a collision with a phonon. The same arguments as we used in the previous section lead to a system very similar to (11-14), except that equation ( 13) now contains a collision term in the right-hand side since gas atom-phonon collisions can occur:

v x ∂ y F + v z ∂ z F - 1 2 ∂ y V # (y, z)∂ vx F - 1 2 ∂ z V # (y, z)∂ vz F = 1 λ * Q m,p [F ], 0 < z. ( 19
)
Again the left-hand side of ( 19) is nothing but the derivative of F along the trajectories of gas atoms, and ( 19) can be written as

d ds F (t, x, y(s), z(s), v(s)) = 1 λ * Q + m,p [F (t, x, y(s), z(s), v(s))] - F (t, x, y(s), z(s), v(s)) τ m,s (E) . ( 20
)
But because of the gain term Q + m,p [F ], it is not possible to get an explicit form for the solution of this equation. However, if we replace this term by Q + m,p [α(t, x)M] where α(t, x) is a constant (with respect to the variables y and z) to be defined later, it is easy to obtain the explicit solution (since Q + m,p [α(t, x)M] and τ ms (E) are constant along the characteristic):

F (t, x, y, 0, v) = exp(-r(y, v))g(t, x, 0, -Λ 2 (y, -v)) + (1 -exp(r(y, v))) α(t, x) exp(-V m )M w (v), (21)
where M w (v) = exp(-m|v| 2 /(2kT w ) and r(y, v) = τ f l (y, -v)/λ * τ ms (E(y, 0, -v)).

With such a choice of the gain term we relax the (local) property of mass conservation of the collision term but we choose α so that the mass flux through the boundary vanishes (which is a global property of mass conservation) (see details in [START_REF] Brull | Gas-surface interaction and boundary conditions for the Boltzman equation[END_REF] ). Finally, after some algebra we obtain the following boundary equation

g(t, x, 0, v)|v z | |vz<0 = v z >0 R(v → v)|v z |g(t, x, 0, v )dv + a(v)σ[g(t, x, 0, v ) |v z >0 ]|v z |M w (v), (22) with R(v → v) = 1 0 exp(-r(y, v)) |v z | |v z | δ(v + Λ 2 (y, -v)) dy, a(v) = 1 - 1 0 exp(-r(y, v)) dy, σ[g(t, x, 0, v ) |v z >0 ] = v z >0 v z 1 - 1 0 exp(-r(y , -v ))dy g(t, x, 0, v )dv -vz<0 v z 1 - 1 0 exp(-r(y, v))dy M w (v)dv .
Note that the above coefficients have the following physical interpretations. First, exp(-r(y, v)) is the fraction of incident gas atoms re-emitted with velocity v which impinged the wall with velocity v = -Λ 2 (y, -v) and that had no collision with phonons, and moreover this term does not depend on the approximation of the gain term. Then 1 0 exp(-r(y, v)) dy = 1a(v) is the fraction of all the incident gas atoms re-emitted with velocity v (whatever incident velocity) and that had no collision with phonons. Finally, a(v) is the fraction of incident gas atoms re-emitted with velocity v and that had a collision with phonons. Moreover, note that for a flat wall (i.e. if we assume that V(x, z) = V(z)), relation (22) yields

g(t, x, 0, v)|v z | |vz<0 = (1 -a(v))|v z |g(t, x, 0, v) + a(v)σ(t, x)|v z |M w (v),
where a(v) is the fraction of incident gas atoms that have been thermalized . This relation is very similar to the classical Maxwell boundary condition but with a coefficient a which depends-as in (22)-explicitly on the velocity v but also on the temperature T w (through τ ms , see (7)) and on the wall morphology (through τ f l ).

The last regime to be studied is when we assume that τ * ms τ * f l , that is to say that the relaxation time is much smaller than the time of flight which implies that all gas atoms entering the surface layer have a collision with a phonon and are thermalized.

Then we get as it is expected the usual boundary condition of diffuse reflection (or total accommodation). This can be obtained either by an asymptotic analysis as in the previous regime studied above or by noting that lim

τms/τ f l →0 a(v) = 1. ( 23
)

III. IMPLEMENTATION OF MLBC IN A DVM CONTEXT

In order to use the boundary conditions proposed in the previous sections, the potential energy surface V # has to be specified. Now we consider that the wall is a solid made up of atoms organized in a periodic structure which is assumed to be defect free. The potential V # is generated by the interaction potentials of these atoms (see [START_REF] Brull | Boundary conditions for the Boltzmann equation for rough walls[END_REF] ), i.e. V # (y, z) = i V i (y, z) where each potential V i is assumed to be a Lennard-Jones potential, with mixed-interaction parameters (in dimensionless form) σ and ˜ , centered at the position of an atom of the lattice

V i (y, z) = V (r i ) = 4˜ ( σ ri ) 12 -( σ ri ) 6 ,
where ri = (yy i ) 2 + (zzi ) 2 , ˜ = 2 /kT and σ = σ/L * . For the computations we will take only into account the atoms located in the periodic cell or in neighboring cells (see Fig. 4). With this definition of the potential energy surface, condition (2), which implies that the wall is non porous, is not satisfied, but this is not a problem for practical applications where the energy of incident gas atoms is not large enough to allow them to penetrate the wall. Note that in this configuration the roughness of the wall is of the order of the lattice period , i.e. at a very small scale (typically a nanometer or less). It is in the order of a standard experimental roughness for silicon channels or tubes of fused silica 12 15 16 .

We first consider the simplest situation of elastic scattering where no atom is thermalized during the gas-surface interaction (see section II B). Then, as we indicated above, for a flat wall the boundary condition is the usual specular reflection. For a periodic rough wall the boundary condition introduced in [START_REF] Brull | Gas-surface interaction and boundary conditions for the Boltzman equation[END_REF] is given by (16-17).

To make the boundary condition ( 16) explicit, we have to compute R(v → v), or 1 0 δ(v + Λ 2 (y, -v)) dy ( for v z > 0 and v z < 0) which is the probability that an impinging gas atom with velocity v is re-emitted with velocity v. At the same time, for a given v , we obtain the scattering pattern of a mono-energetic beam of incident gas atoms by the wall and we can observe the effect of the roughness induced by our model. In the context of a numerical approximation of the Boltzmann (or BGK) equation by a discrete velocity method (DVM) 17 18 we have to compute 1 0 δ(v p + Λ 2 (y, -v q )) dy, for each pair (v p , v q ) in the set of discrete velocities. Since we are in a regime where the gas atom-phonon collisions are neglected (structure scattering), R(v → v) is obtained by computing the trajectories of incident gas

ẏ = v x , vx = -∂ y V(y, z), (24) ż
= v z , vz = -∂ z V(y, z). (25)
It is a Hamiltonian system in which the total energy E(y, z, v) = 1 2 (v 2 x + v 2 z ) + V(y, z) remains constant along the trajectories. It is numerically solved by using the velocity Verlet method which is a classical symplectic integrator. As a consequence of the energy conservation, the norm of the velocity of a gas atom leaving the surface layer is the same as when the gas atom entered the surface layer (since V(y, z = 0) is independent of y), but of course the re-emission angle has no simple relation with the incident angle. More precisely, if the velocity space is discretized by the following grid in polar coordinates:

V = {v k,i = r k (cos θ i , sin θ i ), 1 ≤ k ≤ k max , 1 ≤ i ≤ i max }, then R(v k ,i → v k,i ) = 0 if k = k . On this grid, the distribution function is approx- imated by f k,i ≈ f (v k,i ).
tor is [-θ i -π/(2i max ), -θ i + π/(2i max )],
which is centered around -θ i (where θ i is the ith discrete incident angle used to define the discrete velocity set). Boundary condition ( 16) is then discretized by

| sin θ i |f k,i = imax i =1 | sin θ i |R(v k ,i → v k,i )f k,i r k ∆r∆θ,
where ∆r = r k+1r k and ∆θ = π/i max are the (constant) steps of the velocity grid in the r and θ directions. Note that in this relation, we used the fact that the terms in which k = k are zero.

Note that for a practical use of this boundary condition, for a same wall material, the scattering kernel R(v k ,i → v k,i ) can be computed and stored once for all.

IV. RESULTS AND DISCUSSION

A. Scattering pattern for elastic scattering on a rough wall

As indicated above, when we compute R(v → v) we obtain at the same time the scattering pattern of a mono-energetic beam of incident gas atoms by the wall and we can observe the effect of the roughness induced by our model. In our simulations the atoms of the wall are located in the periodic cell at points (y = 0, z = 2), (y = 1, z = 2). The interaction potential of the wall is made up of the interaction potentials of these two atoms and of the adjacent atoms in the lattice (see Fig. 4). The temperature T w was set to 296 K, L * to 0.4 nm, the molecular mass m = 123 amu 2.04 × 10 -25 kg, so that the reference velocity is v * = 200 m/s and the dimensionless parameters of the LJ potential are set to σ = 1, ˜ = 1 (corresponding to σ = 0.4 nm and = 2.04 10 -21 J ), and L = 2 σ = 0.8 nm.

In the following, we plot for each ith angular sector the probability p i that an incident atom with velocity v inc is re-emitted in this sector. It is computed as follows:

we count the number n i of gas atoms re-emitted in the ith angular sector, and then we set p i = n i /m max , where m max is the number of trajectories. The results presented below are computed with m max = 8001 trajectories, since the results were not significantly modified when we used more trajectories. All trajectories were computed using a velocity Verlet method with a time step ∆t = 2 × 10 -2 fs. For the visualization we take i max = 41 so that each angular sector has an angle of 4.30

•
and there is an average of 195.15 re-emitted atoms per angular sector.

Influence of the angle of incidence

In this simulation we take |v inc | = 300 m/s, and θ inc takes the values (in degrees)

Influence of the energy (or mass) of the incident gas atoms

The trajectory of an incident gas atom through the surface layer depends on the position (y, 0) where it enters the surface layer and on its incident kinetic energy (since V(y, 0) is independent of y). So changing the incident kinetic energy can be interpreted as keeping the same incident velocity but changing the mass of the gas atoms (i.e. changing the gas). The results given below are obtained by setting the angle of incidence to 45 amu or m = 2 amu), the largest peak is obtained for forward scattered gas atoms with a reflection angle very close to the specular reflection angle θ = -45 • . Finally we note that the scattering patterns depend significantly on the mass of the gas atom. This should be taken into account for simulating flows of mixtures of gas with different molecular masses, for instance when considering separation processes through micro-porous materials. Let us recall that the influence of the mass of the gas atom on the gas-surface interaction was shown by experiments in previous papers 12 15 . More precisely experiments suggest that the wall/gas interaction tends to be more diffuse with decreasing molecule mass and that the wall asperities can manage and scatter more easily the light incoming molecules. We note that the results of Fig. 6 are quite in agreement with those remarks .

Influence of of the interaction potential

The role of the parameters of the Lennard-Jones interaction potential is investigaded in Fig. 8. We compare the scattering pattern for an incident beam with θ inc = 45 On Fig. 8, the projections of the polar diagram of the scattering pattern on the plane (v y 1 , v z ) and on the plane (v y 2 , v z ) are shown for the two simulations with θ inc = 90 • and θ inc = 60 • . For θ inc = 90 • , the scattering plane is aligned with the lattice while it is not aligned with the lattice for θ inc = 60 • . The influence of the azimuthal angle is noticeable. As expected the projection on the (v y 1 , v z ) plane is symmetric when the scattering plane is aligned with respect to the symmetry axis of the surface lattice, while it is not at all symmetric when the scattering plane is misaligned (d). Moreover the projections on the (v y 2 , v z ) plane are quite different, showing a large peak of backward scattered gas atoms (around -98 • ) for θ inc = 60 • .

B. Non-elastic scattering on a rough wall

As shown in II C, the fraction of gas atoms that are thermalized during the interaction with the surface is given in our model by

a(v) = 1 - 1 0 exp(-τ f l (y, -v)/(λ * τ ms (E(y, 0, -v))) dy,
and so depends on v, ˜ , and σ, through τ f l , and also on T w , the temperature of the wall, through τ ms . The time of flight τ f l is easily obtained when computing the trajectory of the gas molecules in the surface layer by solving (24-25). The time τ ms is a decreasing function of the wall temperature T w , but is unfortunately more difficult to evaluate. Until now we lack experiments, but approximations could be obtained through molecular dynamics simulations. This is a difficult problem by itself and will be treated in forthcoming publications. Here, we follow a simplified approach 9 that consists in approximating a by

a(v) = 1 - 1 0 exp(-τ f l (y, -v)/(λ * τ mean ms (T w )) dy, (26)
where τ mean ms (T w ) is an average relaxation time due to gas atom-phonon collisions. In order to investigate the effects of the temperature of the solid, we simulate a(v) as a function of τ mean ms . Since τ mean ms ≥ 10 -13 s (the characteristic period of the thermal vibration of atoms in a solid 9 ), and takes values roughly in the range 10 -13 -10 -12 s (see 9 11 20 ), we studied the variation of a(v) for 0.25 10 -12 s ≤ τ mean ms ≤ 1.4 10 -11 s (Fig. 10a), for a beam of impinging gas atoms with m = 123 amu, θ inc = 45 • and λ * = 2. We use a simple numerical approximation of the integral of (26) to get:

a(v) ≈ 1 - mmax m=1 exp(-τ f l (y m , -v)/(λ * τ mean ms (T w )) ∆y,
where τ f l (y m , -v) is the time for a particle that enters into the surface layer in y m = m∆y with a velocity -v to leave this surface: it was computed exactly when trajectories equations (24-25) were solved in section IV A.

As expected, the fraction of incident gas atoms that have been thermalized decreases. For τ mean ms = 0.25 10 -12 s, the fraction is very closed to 1, which means that we are in the regime of total accommodation, while for τ mean ms = 1.4 10 -11 s this fraction is very small so that the scattering pattern is mainly governed by the morphology of the wall ("structure scattering"). Figure 10b Finally the influence of the incident angle on the coefficient a is shown in Table I where the simulations were performed with λ * = 2 and τ mean ms = 0.525 10 -12 s . It is clear that the smaller the incident angle, the larger the fraction of thermalized gas atoms, as a consequence of larger times of flight.

V. CONCLUSION

We have proposed new boundary conditions allowing to take into account at the kinetic level the effect of the van der Waals forces and of the periodic roughness of the wall at the atomic scale. This model relies on a two scale theoretical approach whose mathematical framework has been introduced in a previous paper. We have proposed a way to include some physics of gas-surface interaction by computing the trajectories of the impinging gas atoms through a periodic cell of the surface layer in which gas molecules interact with surface atoms by a simplified PES. Then it is possible to get an approximation of the scattering kernel of the new boundary condition suitable for a DVM numerical approach of the Boltzmann equation. We have shown the scattering pattern of a mono-energetic beam of impinging gas atoms reflected by a wall with nanoscale roughness. Several simulations on a simplified 2D configuration or on a 3D configuration prove that the model brings at the kinetic scale information on the roughness of the wall at the atomic scale which are not present in usual boundary conditions. 2D simulations of elastic scattering show that the scattering pattern is strongly dependent on several parameters such as the incident angle of the impinging gas atoms, their energy or the mixed-interaction parameters of the Lennard-Jones potential. Moreover 3D computations prove the influence of the azimuthal angle of the scattering plane. For non-elastic scattering a simplified form of the new boundary condition shows that the temperature of the wall and the incident angle of the impinging gas atoms have a noticeable influence on the fraction of thermalized gas atoms. But further investigations are necessary to obtain a more precise evaluation of the relaxation time of gas atoms by the phonons. This needs some more experiments and/or numerical Molecular Dynamics simulations and will be treated in forthcoming publications. All those simulations are promising and prove that our approach is simple enough to be numerically tractable, while it is able to include enough physics of gas-surface interaction so as to improve kinetic simulations of micro or nano-flows in narrow channels, and hence could be a useful tool in this context.

FIG. 3 .

 3 FIG.3. One trajectory (and its corresponding reversed trajectory in blue) and the mapping Λ.

FIG. 4 .

 4 FIG. 4. Periodic cell of the surface layer (2D case). The interaction potential of the colored atoms are taken into account

  For a given incident velocity v inc = r inc (cos θ inc , sin θ inc ), several trajectories are computed for various equi-distributed initial condition s (y m , v inc y , v inc z ), with y m = m∆y, 0 ≤ m ≤ m max , ∆y = 1/m max . For each trajectory we obtain the angle of reflection θ ref of the gas atom in the flow (v ref = r inc (cos θ ref , sin θ ref )). The interval [-180 • , 0] of reflected angles is divided into i max /2 angular sectors (like the interval of incident angles). The ith angular sec-

22. 5 •

 5 , 45 • , 67.5 • , 90 • .

Figure 5 FIG. 5 .

 55 Figure 5 shows very clearly that the scattering pattern is quite different from the usual specular reflection for any incident angle. For θ inc = 90 • the diagram is symmetric and shows several peaks with a concentration of re-remitted particles near the specular reflection angle (90 • ). For the other values of θ inc , the diagram shows a lobal distribution of backward scattered gas atoms with a peak more or less close to the specular reflection angle and a high but narrow peak of backward scattered gas atoms. More precisely, for θ inc = 22.5 • , the lobe shows several local peaks between -45 • and -75 • (for a specular angle of reflection equal to -22, 5 • ), and the peak of backward scattered gas atoms is around -107 • . For θ inc = 45 • , the lobe shows two local peaks around -40 • and -50 • for a specular reflection angle equal to -45 • and the peak of backward scattered gas atoms is between -121 • and -125 • . For θ inc = 67.5 • , the lobe shows a peak at -70 • , for a specular reflection angle equal to 67.5 • and the peak of backward scattered gas atoms is at -147 • .

  • and setting the norm of the velocity to |v inc | = 6v * = 1200 m/s. Figure 6 shows the computed polar diagrams of the scattering patterns for gas atoms with masses 123, 60, 16 and 2 amu. For heavy atoms (m = 123 amu) we observe a small peak in the polar diagram for θ = -20 • and a larger one for θ = -85.6 • . For m = 60 amu the larger peak is for θ = -90 • and a smaller one for θ = -37.3 • . The polar diagram corresponding to m = 16 amu shows a backward scattering with a large peak at θ = -103.1 • and a forward scattering with a smaller peak at θ = -46 • . Finally for the lightest atoms (m = 2 amu) the diagrams shows also a backward scattering with a small peak at θ = -103.1 • and a forward scattering with larger peak at θ = -46 • . Note that for light atoms (m = 16

,

  |v inc | = 1200 m/s with an interaction potential V 1 of parameters σ = 0.4 nm and = 4.08 10 -21 J and for an interaction potential V 2 with parameters σ = 0.6 nm and = 8.16 10 -22 J (see Fig.7 (a)).Both polar diagrams show two peaks but at very different angles: -20 • and -85 • for the first and -35 • and -55 • (closer to the specular reflection angle) for the other one. Those results indicates that the scattering pattern also strongly depends on the material of the wall.

FIG. 6 .FIG. 7 .

 67 FIG. 6. Polar diagrams of the scattering pattern for impinging gas atoms with an incident angle 45 • , a velocity |v inc | = 1200 m/s and for gas masses 123 (a), 60 (b), 16 (c), and 2 amu (d).

FIG. 8 .

 8 FIG. 8. Polar diagrams of the scattering pattern for impinging molecules with an incident angle 45 • , a velocity |v inc | = 1200 m/s, for a potential V 1 with σ = 0.4 nm and = 2.04 10 -21 J (a) and for a potential V 2 with σ = 0.6 nm and = 1.63 10 -21 J (b).

FIG. 9 .

 9 FIG. 9. Projection of the polar diagram of the scattering pattern on (v y 1 , v z ) for θ = 90 • (a), on (v y 2 , v z ) for θ = 90 • (b), on (v y 1 , v z ) for θ = 60 • (c), and on (v y 2 , v z ) for θ = 60 • (d).

  shows the dependence of the fraction a(v) on the the velocity. We simulate |v inc | from 200 m/s up to 10000 m/s for τ mean ms = 0.525 10 -12 s and m = 123 amu. It appears that for small velocities, the coefficient a ≈ 0.8 while for |v inc | = 2400 m/s a ≈ 0.13, and we are close to an elastic scattering, which is a consequence of the smaller times of flight τ f l (y) of each trajectory.

  FIG. 10. Fraction a(v) of thermalized molecules, for (a) |v inc | = 1200 m/s, θ inc = 45 • , λ * = 2, as a function of τ mean ms , and for (b) θ inc = 45 • , τ mean ms = 5.25 10 -12 s, λ * = 2, as a

TABLE I .

 I Fraction of thermalized gas atoms as a function of the incident angle incident angle. fraction of thermalized gas atoms.

	22.5 •	0.3736
	45 •	0.23
	67.5 •	0.1999
	90 •	0.1906


Jessy Mallet
email: mallet@celia.u-bordeaux1.fr

Stéphane Brull
email: stephane.brull@math.bordeaux1.fr

Bruno Dubroca
email: dubroca@celia.u-bordeaux1.fr

GENERAL MOMENT SYSTEM FOR PLASMA PHYSICS BASED ON MINIMUM ENTROPY PRINCIPLE

Keywords: 2000 Mathematics Subject Classification. 35Q84, 35Q83, 65Z05 entropy minimization, Landau-Fokker-Planck equation, N moment system, entropic scheme, moment closure, entropic average

In plasma physics domain, the electrons transport can be described from kinetic and hydrodynamical models. Both methods present disadvantages and thus cannot be considered in practical computations for Inertial Confinement Fusion (ICF). That is why we propose in this paper a new model which is intermediate between these two descriptions. More precisely, the derivation of such models is based on an angular closure in the phase space and retains only the energy of particles as a kinetic variable. The closure of the moment system is obtained from a minimum entropy principle. The resulting continuous model is proved to satisfy fundamental properties. Moreover the model is discretized w.r.t the energy variable and the semi-discretized scheme is shown to satisfy conservation properties and entropy decay.

Introduction.

A variety of classical problems in kinetic theory leads to use the Fokker-Planck-Landau equation to describe the evolution of different species of charged particles undergoing binary collisions ( [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF][START_REF] Delcroix | Physique des Plasmas[END_REF]). For these collisions, the interaction potential is the long-range Coulomb interactions. More precisely the solutions of the kinetic equations are non-negative distribution functions f α (t, x, v) specifying the density of each species α with velocity v at time t and position x.

Here the plasma consists of electrons and ion species (see [START_REF] Tzoufras | A Vlasov-Fokker-Planck code for high energy density physics[END_REF]). In this paper, the plasma is studied during the time scale corresponding to the electron frequency. This time scale being small compared to the characteristic time of ions motion, the ions can be considered as immobile. To approximate the solution of such problems, many computational methods have been developed up to now.

In [START_REF] Buet | Numerical analysis of conservative and entropy schemes for the Fokker-Planck-Landau equation[END_REF][START_REF] Buet | Fast algorithms fot numerical, conservative, and entropy approximations of the Fokker-Planck-Landau equation[END_REF][START_REF] Crouseilles | Numerical approximation of collisional plasmas by high order methods[END_REF], deterministic schemes have been developed to approximate the solution of the Fokker-Planck-Landau equation by using a phase space grid. Then in [START_REF] Duclous | High order resolution of the Maxwell-Fokker-Planck-Landau model intended for ICF application[END_REF], the evolution of electromagnetic fields was inserted in the scheme to describe exactly the electrons transport. This scheme presents many fundamental properties such as conservation laws and positivity of the distribution function. However kinetic models for electrons transport are too much expensive to be used for Inertial Confinement Fusion (ICF). To reduce the computational time, plasmas can be described by fluid models. For example in [START_REF] Crispel | A plasma expansion model based on the full Euler-Poisson system[END_REF][START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF][START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasi-neutral limit[END_REF], the authors consider a bi-fluid compressible Euler model coupled with the Poisson equation. Moreover for ITER application, the isothermal two-fluid Euler-Lorentz system coupled with a quasi-neutrality constraint has been studied in [START_REF] Brull | An asymptotic preserving scheme for a bi-fluid Euler-Lorentz system[END_REF] by introducing an asymptotic preserving scheme. However, for the new high energy target drivers, the kinetic effects are too important to neglect them.

In the present paper, we propose an alternative approach by considering an intermediate description between fluid and kinetic level like in [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Maxwell-Fokker-Planck equation of electrons[END_REF]. The velocity variable is written in spherical coordinates and the model is written by considering moment systems with respect to the angular variable. But the electron/electron collision operator being non linear, the moments extraction is complicated. That is why physicists approximate this operator by assuming that the main contribution of the distribution function comes from its isotropic part ( [START_REF] Yu | Conservative finite-difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF]). But this approached model does not conserve the realizability domain defined as the set of vectors that are the moments of positive distribution functions.

There exists many moment models whose difference comes from the choice of the closure which approximates the full distribution function. This closure is essential to get a closed moment system and to assure some physics properties. Indeed, a well-known moment model called P N (briefly reviewed in [START_REF] Case | Linear transport theory[END_REF]) leads to solutions that are inconsistent with a positive concentration of particles. For example, the P N model used in [START_REF] Mcclarren | On solutions to the Pn equations for thermal radiative transfert[END_REF] does not satisfy positivity of the underlying distribution function of electrons and entropy dissipation. This comes from the definition of the closure that is based on a truncation of the spherical harmonic expansion w.r.t the angular variable. For the P N model, modifications are proposed in ( [START_REF] Hauck | Positive P N closure[END_REF]) to obtain a non-negative distribution function. For example, in [START_REF] Larsen | The P N theory as an asymptotic limit of transport theory in planar geometry -I: Analysis[END_REF][START_REF] Larsen | The P N theory as an asymptotic limit of transport theory in planar geometry -II: Numerical results[END_REF], the authors propose an asymptotic limit approximation to prescribe boundary conditions for the neutronic transport equation.

In this context, the motivation of our current work is to use a moments closure based on an entropy minimization principle which inherits many of the fundamental properties such as the conservation laws and entropy dissipation. This closure leads to a new moment model called M N model. This model has been firstly introduced in [START_REF] Dubroca | Entropic moment closure hierarchy for the radiative transfert equation[END_REF] in a reduced form to find an approximation of radiative transfer and energy evolution equation. The subsequent system of conservation laws is closed by minimizing the radiative entropy (see also [START_REF] Frank | Partial moment entropy approximation to radiative transfer[END_REF] for extension). Moreover, a moment model M 1 has been proposed in [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Maxwell-Fokker-Planck equation of electrons[END_REF] by considering the three first moments. The benefit of the M 1 model lies in that it satisfies fundamental properties. Indeed mass and energy are conserved, the distribution function stays positive and the entropy decreases. Nevertheless proofs presented in [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Maxwell-Fokker-Planck equation of electrons[END_REF] are very specific to the M 1 model. Indeed theorems and proofs are based on the knowledge of known coefficients to define the distribution function f . Here these coefficients are unknown which implies to use a different approach.

In kinetic theory, the idea of using minimization problem under moments constraints has been analyzed by Levermore in the framework of gas dynamics ( [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF]). The main motivation was to construct BGK models ( [START_REF] Bathnagar | A model for collision processes in gases[END_REF]) leading to the correct Prantl number at the Navier-Stokes level. Moreover the well-posedness of such problems has been clarified later by Junk and Schneider ([25,[START_REF] Schneider | Entropic approximation in kinetic theory[END_REF]). In the present case, the aim is different. Here the energy of particles is a free parameter. Then we integrate only the kinetic equation expressed in spherical coordinate with respect to the angle variable and we return only the energy of particles as a kinetic variable in the frame of ions considered completely frozen.

At the numerical point of view, the entropic average used to define the distribution function on interfaces allows to show that the semi-discretized scheme is entropic.

The remainder of the paper is organized as follows. In section 2, we briefly review some backgrounds about the collision operator properties and present the kinetic equation considered. In section 4, we introduce a new electron/electron collision operator that is devoted to approach the full Landau collision operator. This new model is based on a linearization of the Landau operator around the equilibrium state of the electron/ion collision operator. Next, this new operator is proved to preserve the realizability domain. Section 5 deals with the derivation of a semi-discrete scheme, where the time is kept as a continuous variable whereas the energy variable is discrete. The definition of the approximated solution on the dual mesh through an entropic average allows to get the entropic dissipation property of collision operator. Next, the N-moment model is derived from the continuous kinetic equation in section 6 and is shown to preserve previous property. In section 7, the semi-discrete scheme for the moment system is proved to be entropic. Moreover we present at the end of this part by a test case, showing the interest of the M 2 and M 3 models compared to the M 1 model. Comparisons with P N models are also given. Finally we finish in the last section by conclusions to this work.

2.

Presentation of the kinetic model. In this section we recall some backgrounds about classical kinetic models for plasma physics where the ions are considered as frozen. Firstly, the general kinetic model is presented. Then some fundamental properties of collision operators are explained.

2.1. Classical kinetic models. In this paper, the time evolution of the electrons is describded by a distribution function f (t, x, v), x ∈ R 3 , v ∈ R 3 solution of a kinetic equation.

hydrodynamic quantities.

From the distribution function, we define the density, the macroscopic velocity and the temperature of the electrons as follows

n = R 3 f (t, x, v)dv, nu = R 3 vf (t, x, v)dv, T = 2 3 R 3 (v -u) 2 f (t, x, v)dv.
2.1.2. Classical kinetic equations. In the present paper, the distribution function f satisfies the kinetic equation

∂ t f + v∂ x f = C(f, f ) , (2.1)
where the collision operator C(f, f ) introduced in [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF] is defined by

C(f, f ) = C ee (f, f ) + C ei (f ) .
The operator C ee stands for the electron-electron collision operator

C ee (f, f ) = ∇ v . R 3 Φ(V )[f (v ′ )∇ v f (v) -f (v)∇ v ′ f (v ′ )]dv ′ , (2.2)
where V = v -v ′ is the relative velocity of electrons and Φ is an operator acting on the relative velocity V

Φ(V ) = 1 |V | 3 (|V | 2 Id -V ⊗ V ) , (2.3)
where Id is the unit tensor.

The electron-ion collision operator C ei (f ) is defined by

C ei (f ) = ∇ v • [Φ(V )∇ v f (v)] . ( 2
R 3 C ee (f, f )   1 v |v| 2   dv = 0, t ≥ 0 ,
while the electron-ion collision operator satisfies only mass and energy conservation

R 3 C ei (f ) 1 |v| 2 dv = 0, t ≥ 0 .
They both dissipate the entropy i.e. C ei (f ) log f dv ≤ 0 and

R 3 C ee (f, f ) log f dv ≤ 0 , (2.5)
which implies that the Boltzmann entropy

H(f ) = R 3 (f log f -f ) dv is a Lyapounov function for (2.1). Property 1. [(i)] 1.
The equilibrium state of the electron-ion collision operator C ei (i.e.

C ei (f ) = 0) is given by the set of isotropic functions f = f (|v|) . 2. The equilibrium state of the electron-electron collision operator C ee (i.e. C ee (f, f ) = 0) is given by the Maxwellian distribution function such that

∃(n, m e , T, u) ∈ R 3 + × R 3 / f (v) = n m e 2πk B T 3 2 exp - m e (v -u) 2 2k B T ,
where k B is the Boltzmann constant, n is the density, T is the temperature and u e represents the mean velocity. 3. The equilibrium state considered by both collision operators is given by the isotropic Maxwellian distribution function such that

∃(m e , n, T ) ∈ R 3 + /f (v) = n m e 2πk B T
f i = 2πζ 2 1 -1 f (ζ, µ)µ i dµ = ζ 2 f µ i ∀i ∈ {0, N }, (3.1)
where . is defined for any function Ψ by

Ψ = 2π 1 -1 Ψ(µ)dµ .
We use also the following notation

F i (ζ) = f i (ζ) ζ 2 . (3.2)
This means that F 0 represents the isotropic part of the distribution function f (F 0 = f ). In this paper, moment systems will be constructed from kinetic equations. One fundamental property for these systems is the preservation of the realizability domain

A =      g =    g 0 . . . g N    ∈ R N / ∃ g ≥ 0 ∈ L 1 ([-1, 1]) and g i = ζ 2 µ i g(µ) ∀i ∈ {0, N }      . (3.3)
3.3 garantees that the moments can be recovered from a nonnegative distribution function.

3.2. The M N closure. One important step of the present work, is the construction of a moment system from a kinetic equation. But the extraction of these moments leads to system of equations having more unknowns than equations because a supplementary moment appears. That is why, a closure must be defined such that the highest moment writes in function of the previous ones. The closure chosen is the M N closure that is based on an entropy minimization principle. The advantage of such a closure compared to the classical P N closure is the nonnegativity of the underlying distribution function. More precisely, the M N closure is obtained by solving the following entropy minimization problem:

f i (ζ), i ∈ {0, N } being given min g≥0 {H(g), ζ ∈ R + , / 1 -1 µ i g(µ, ζ)dµ = f i (ζ), i ∈ {0; N }}. (3.4)
Therefore according to [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF][START_REF] Minerbo | Maximum entropy Eddington Factors[END_REF], if

(f i (ζ)) i∈{0;N } ∈ A, then f solution to (3.4) reads f (t, x, ζ, µ) = exp (α(t, x, ζ) µ) , (3.5
)

with α =    α 0 . . . α N    and µ =    µ 0 . . . µ N   .
3.3. Conservation of the realizability domain. Consider here the case of a space homogeneous situation with one species of particle, i.e. the plasma consists only of electrons. In this case, equation (2.1) is reduced to

∂ t f = C ee (f, f ) . (3.6)
The moments extraction for the electron-electron collision operator C ee (f, f ) is very complicated because of its nonlinearity. That is why in Plasma physics, classical approximations for the operator lead to consider that the main contribution for the electron-electron collision operator comes from the isotropic part of the distribution function. This means that the collision operator C ee (f, f ) is approached by

Q 0 ee = 1 ζ 2 ∂ ζ ζ ∞ 0 J(ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ (F 0 (ζ)) -F 0 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ , (3.7) with
J(ζ, ζ ′ ) = 2 3 inf 1 ζ 3 , 1 ζ ′3 ζ ′2 ζ 2 . (3.8)
This approximation (3.7) denoted by Q 0 ee has been presented in [START_REF] Yu | Conservative finite-difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF] and used in [START_REF] Buet | Conservative and Entropy Decaying Numerical Scheme for the isotropic Fokker-Planck-Landau Equation[END_REF][START_REF] Dellacherie | Contribution à l'analyse et à la simulation numérique des équations cinétiques décrivant un plasma chaud[END_REF][START_REF] Dellacherie | Sur un schéma numerique semi-discret appliqué à un opérateur de Fokker-Planck isotrope[END_REF] to study the homogeneous Fokker-Planck-Landau equation for isotropic distribution functions. However, the following example 1 shows that this model does not preserve the realizability domain A.

Indeed, consider the distribution function as in (3.5) and extract the two first moments of (3.6) with respect to µ, where C ee (f, f ) has been replaced by Q 0 ee . Hence, we get

   ∂ t f 0 = Q 0 ee , ∂ t f 1 = 0 . (3.9)
For more details, the expression (3.7) and the derivation of the system (3.9) are explained in Appendix A.

In the particular case of the N = 1, the realizability domain A is shown to be equal to

B = g = g 0 g 1 ∈ R 2 , g 0 > 0 and |g 1 | < g 0 ∪ {(0, 0)} .
Example 1. Let us chose the following initial data for f 0 and f 1

f 0 (t = 0) = 1 3 χ [0,3] (ζ) and f 1 (t = 0) = 1 4 χ [0,3] (ζ) .
Because of the electron-electron collision operator effect, f 0 converges to a centered Maxwellian. Moreover, the M 1 model preserves mass and energy. So the obtained Maxwellian presents the same mass ( ∞ 0 f 0 dζ = 1) and the same energy

( ∞ 0 f 0 ζ 2 dζ = 3
) as the initial distribution function. So the asymptotic steady state for f 0 reads

f 0 = 2 π exp -ζ 2 2 ζ 2 .
Therefore figure 1 shows that the realizability domain is not preserved anymore. 4. New Kinetic continuous model.

4.1.

Approximation of the collision operator C ee . In order to preserve the realizability domain (proof in section 6.2), we consider a new collision operator based on a linearisation of C ee around the equilibrium state of C ei . The expression of this new collision operator is given by equation (4.1).

Approximation of C ee The electron-electron collision operator Q ee (f ) and the electron-ion collision operator Q ei (f ) are given by

Q ee (f ) = 1 ζ 2 ∂ ζ ζ ∞ 0 J(ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ f (ζ) -f (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ , (4.1) where J(ζ, ζ ′ ) is given in (3.8).
For C ei , we do not use any approximation, since this operator is already linear. To be consistent with the notation Q ee (f ), we use the notation Q ei instead of C ei . Its expression can be rewritten

Q ei (f ) = 1 ζ 3 ∂ ∂µ (1 -µ 2 ) ∂f ∂µ . (4.2)
The computation of (4.2) in spherical coordinate is established in Appendix B. So equation (2.1) can be approached in spherical coordinates by

∂ t f + ζµ∂ x f = Q(f ) , (4.3)
where

Q(f ) = Q ee (f ) + Q ei (f ) and Q ee , Q ei are defined by (4.1, 4.2). Remark 1. If f is isotropic, i.e. f = F 0 , we obtain Q(f ) = Q(F 0 ) = 1 -1 C ee (F 0 , F 0 )dµ = 1 -1 C ee (f, f )dµ .
Now, let us check that the new model conserves fundamental properties.

4.2.

Properties of the model. In this subsection, we present the fundamental properties satisfied by the new model, that is the conservation laws and a H-theorem.

Proposition 1. The operator Q(f ) satisfies mass and energy conservation properties and the entropy dissipation property, i.e.

∞ 0 ζ 2 Q(f ) 1 ζ 2 dζ = 0 and ∞ 0 ζ 2 Q(f ) log f dζ ≤ 0 .
Proof. Firstly, in order to recover the conservations property, consider the quantity

∞ 0 ζ 2 Q(f ) 1 ζ 2 dζ = ∞ 0 ζ 2 Q ee (f ) 1 ζ 2 dζ + ∞ 0 ζ 2 Q ei (f ) 1 ζ 2 dζ .
(4.4) For the sake of clarity, we introduce the notation

Q 0 ee = ζ 2 Q ee (f ) .
The first term of the right hand side of equation (4.4) can be rewritten as

∞ 0 Q 0 ee 1 ζ 2 dζ = ∞ 0 ∂ ζ ζ ∞ 0 J(ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ F 0 (ζ) -F 0 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ 1 ζ 2 dζ . So, ∞ 0 Q 0 ee dζ = 0. Moreover, since ∞ 0 Q 0 ee ζ 2 dζ = -2 ∞ 0 ∞ 0 ζ 2 ζ ′2 J(ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ F 0 (ζ) dζ ′ dζ + 2 ∞ 0 ∞ 0 ζ 2 ζ ′2 J(ζ, ζ ′ ) F 0 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) dζ ′ dζ , we get from Fubini theorem ∞ 0 Q 0 ee ζ 2 dζ = 0 .
Besides, we have Q ei (f ) = 0. Therefore the conservation properties follow. Now we prove the entropy dissipation property of the whole collision operator Q. By using a Green formula we obtain easily that Q ei dissipates entropy. Moreover

∞ 0 ζ 2 Q ee (f ) log f dζ = - ∞ 0 ∞ 0 ζ ′ 2 ζ J (ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ f (ζ) -f (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) 1 f ∂ ζ f dζdζ ′ . (4.5)
Equation (4.5) can be rewritten in terms of f as

∞ 0 ζ 2 Q ee (f ) log f dζ = - 1 -1 1 -1 ∞ 0 ∞ 0 ζ ′ 2 ζ J(ζ, ζ ′ ) f (ζ ′ ) 1 ζ ∂ ζ f (ζ) -f (ζ) 1 ζ ′ ∂ ζ ′ f (ζ ′ ) 1 f ∂ ζ f dζdζ ′ dµdµ ′ = - 1 -1 1 -1 ∞ 0 ∞ 0 ζ ′ 2 ζ 2 J(ζ, ζ ′ )f (ζ)f (ζ ′ ) 1 ζ ∂ ζ (log f (ζ)) - 1 ζ ′ ∂ ζ ′ (log f (ζ ′ )) 1 ζ ∂ ζ log f (ζ)dζdζ ′ dµdµ ′ . Moreover the changing of variables (ζ, ζ ′ ) → (ζ ′ , ζ) leads to ∞ 0 ζ 2 Q ee (f ) log f dζ = - 1 2 1 -1 1 -1 ∞ 0 ∞ 0 ζ ′ 2 ζ 2 J(ζ, ζ ′ )f (ζ)f (ζ ′ ) 1 ζ ∂ ζ (log f (ζ)) - 1 ζ ′ ∂ ζ ′ (log f (ζ ′ )) 2 dζdζ ′ dµdµ ′ , (4.6)
and the entropy dissipation property follows.

Remark 2. The collision operator Q ee does not preserve the impulsion. However, as the whole collision operator C does not preserve also the impulsion, the conservation properties of C are not affected.

5.

Semi-discretized kinetic equation. In this section we propose an energy discretization for the new continuous model defined in section 4 such that the fundamental properties of collision operator are conserved. We firstly discretize the kinetic equation and then take moments in section 7. The key point of the following scheme is the approximation of the distribution function f on the dual mesh, through an entropic average. This average leads to the entropy dissipation property for the discretized collision operator.

Energy discretisation.

Let us define the primal mesh M, for the energy variable ζ, decomposed into a family of rectangles

M j-1 2 =]ζ j-1 , ζ j [, j ∈ {1
, m}, and m ∈ N corresponds to the number of points which discretize the energy domain. ∆ζ j = ζ j -ζ j-1 represents the discretization step, which can be variable.

We denote by D its associated dual mesh consisting of cells

D j =]ζ j-1 2 , ζ j+ 1 2 [, ζ j-1 2 = (j - 1 2 ) ∆ζ j , j ∈ {1; m}, ζ -1 2 = 0. The step ∆ζ j+ 1 2 of D j writes ∆ζ j+ 1 2 = ζ j+ 1 2 -ζ j-1 2
. Let h j be an approximation of h(ζ j ) for all distribution function h and h j+ 1 2 an approximation of h(ζ j+ 1 2 ). The discrete form of (4.3) reads for any j ∈ {1; m},

∂ t f j + ζ j µ ∂ x f j = Q j , Q j = Q ee,j + Q ei,j , (5.1)
where the expressions of Q ee,j and Q ei,j are given in Definition 1.

Definition 1. The collision operators (Q ee,j ) j∈{1;m} and (Q ei,j ) j∈{1;m} are respectively defined by

           Q ee,j = G j+ 1 2 -G j-1 2 ζ 2 j ∆ζ j , Q ei,j = 1 ζ 3 j ∂ ∂µ (1 -µ 2 ) ∂f j ∂µ , (5.2)
with

G j+ 1 2 = ζ j+ 1 2 m k=0 J(ζ j+ 1 2 , ζ k+ 1 2 ) F 0 k+ 1 2 1 ζ j+ 1 2 f j+1 -f j ∆ζ j+ 1 2 - 1 ζ k+ 1 2 f j+ 1 2 F 0 k+1 -F 0 k ∆ζ k+ 1 2 ζ 2 k+ 1 2 ∆ζ k+ 1 2 , (5.3)
where J has been defined in (3.8) and f j+ 1 2 is defined by the following entropic average

     f j+ 1 2 = fj+1-fj log |fj+1|-log |fj | if f j+1 = f j , f j+ 1 2 = f j+1 else .
(5.4)

Remark 3. This entropic average has already been considered for an isotropic distribution function by Dellacherie ([15]) to construct an entropic scheme. Next, in [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Maxwell-Fokker-Planck equation of electrons[END_REF], this approach has been generalized to the M 1 model.

Remark 4. G j+ 1 2 can be simplified as

G j+ 1 2 = A j+ 1 2 f j+1 -f j ∆ζ j+ 1 2 + B j+ 1 2 f j+ 1 2
where

A j+ 1 2 = ζ 2 j+ 1 2 m k=0 min 1 ζ 3 j+ 1 2 , 1 ζ 3 k+ 1 2 F 0 k+ 1 2 ζ 4 k+ 1 2 ∆ζ k+ 1 2 , (5.5)
B j+ 1 2 = -ζ 3 j+ 1 2 m k=0 min 1 ζ 3 j+ 1 2 , 1 ζ 3 k+ 1 2 ζ 3 k+ 1 2 (F 0 k+1 -F 0 k ) . (5.6)
Lemma 1. B j+ 1 2 defined in (5.6) can be simplified into

B j+ 1 2 = ζ 3 j+ 1 2 j k=0 (ζ 3 k+ 1 2 -ζ 3 k-1 2 ) F 0 k .
The property is left to appendix D.

5.2.

Properties of the semi-discretized scheme.

Proposition 2. (Q j ) j∈{1;m} satisfies the following fundamental properties.

1) The operator (Q j ) j∈{1;m} satisfies mass and energy conservations property i.e.

m j=0 ζ 2 j Q j 1 ζ 2 ∆ζ j = 0 .
2) By defining f j+ 1 2 through the entropic average given by (5.4), the operator (Q j ) j∈{1;m} satisfies the entropy dissipation property

m j=0 ζ 2 j Q j log f j ∆ζ j ≤ 0 .
Proof. First we aim to prove the mass conservation. A simple computation gives

m j=0 ζ 2 j Q j ∆ζ j = G m+ 1 2 -G -1 2 .
We suppose G as a compact support function, so

G m+ 1 2 = 0. Besides G -1 2 = 0 because ζ -1 2 = 0. Therefore we get m j=0 ζ 2 j Q j ∆ζ j = 0 .
Then we will show the conservation of energy which can be rewritten into

m j=0 ζ 2 j Q j ζ 2 j ∆ζ j = - m-1 j=0 G j+ 1 2 ζ 2 j+1 -ζ 2 j .
(5.7)

By using the expression of G j+ 1 2 given in (5.3), equation (5.7) can be simplified into

m j=0 ζ 2 j Q j ζ 2 j ∆ζ j = - m-1 j=0 ζ j+ 1 2 m k=0 J(ζ j+ 1 2 , ζ k+ 1 2 ) F 0 k+ 1 2 1 ζ j+ 1 2 f j+1 -f j ∆ζ j+ 1 2 - 1 ζ k+ 1 2 f j+ 1 2 F 0 k+1 -F 0 k ∆ζ k+ 1 2 ζ 2 k+ 1 2 ∆ζ k+ 1 2 ζ 2 j+1 -ζ 2 j . (5.8) Since ζ 2 j+1 -ζ 2 j = 2ζ j+ 1 2 ∆ζ j+ 1 2 , equation (5.8) can be rewritten into m j=0 ζ 2 j Q j ζ 2 j ∆ζ j = - 1 -1 1 -1 m j,k=0 ζ 2 j+ 1 2 ζ 2 k+ 1 2 ∆ζ j+ 1 2 ∆ζ k+ 1 2 J(ζ j+ 1 2 , ζ k+ 1 2 )f k+ 1 2 1 ζ j+ 1 2 f j+1 -f j ∆ζ j+ 1 2 dµdµ ′ + 1 -1 1 -1 m j,k=0 ζ 2 j+ 1 2 ζ 2 k+ 1 2 ∆ζ j+ 1 2 ∆ζ k+ 1 2 J(ζ j+ 1 2 , ζ k+ 1 2 )f j+ 1 2 1 ζ k+ 1 2 f k+1 -f k ∆ζ k+ 1 2 dµdµ ′ ,
and the conservation of energy follows by exchanging indexes k and j. Next, we prove that the numerical operator Q j is entropic. Firstly, by arguing like for the continous case, we get that Q ei,j dissipates the entropy. Moreover

m j=0 ζ 2 j Q ee,j log f j ∆ζ j = - m-1 j=0 log f j+1 f j G j+ 1 2 .
(5.9)

From the definition of G j+ 1 2 , equation (5.9) can be written as

m j=0 ζ 2 j Q ee,j log f j ∆ζ j = - 1 -1 1 -1 m j,k=0 ζ 2 j+ 1 2 ζ 2 k+ 1 2 ∆ζ j+ 1 2 ∆ζ k+ 1 2 J(ζ j+ 1 2 , ζ k+ 1 2 )f k+ 1 2 f j+ 1 2 1 ζ j+ 1 2 f j+1 -f j f j+ 1 2 ∆ζ j+ 1 2 - 1 ζ k+ 1 2 f k+1 -f k f k+ 1 2 ∆ζ k+ 1 2 log f j+1 -log f j ζ j+ 1 2 ∆ζ j+ 1 2 dµdµ ′ .
(5.10) By using the change of variables (µ, µ ′ ) → (µ ′ , µ), the change of index (j, k) → (k, j) and the entropic average (5.4) defining f j+ 1 2 , equation (5.10) can be rewritten as

m j=0 ζ 2 j Q ee,j log f j ∆ζ j = - 1 2 1 -1 1 -1 m j,k=0 ζ 2 j+ 1 2 ζ 2 k+ 1 2 ∆ζ j+ 1 2 ∆ζ k+ 1 2 J(ζ j+ 1 2 , ζ k+ 1 2 ) f k+ 1 2 f j+ 1 2 1 ζ j+ 1 2 log f j+1 -log f j ∆ζ j+ 1 2 - 1 ζ k+ 1 2 log f k+1 -log f k ∆ζ k+ 1 2 2 dµdµ ′ ,
and the entropic dissipation property follows.

6. Continuous moment model. In this section, we explain the construction of the M N model obtained from the kinetic model given in section 4. Besides collision operator properties are shown to be still preserved by the moments extraction.

6.1. N-moment model for the new kinetic equation. Firstly the M N model constructed from the previous kinetic equation (4.3) is presented. This system is established by using a minimum entropy principle for the angle variable, keeping the energy of particles as a kinetic variable and a moments extraction.

The expression Q i ee and Q i ei are computed in the following property. Property 2. The moments of Q ee are given by

Q i ee (f ) = ∂ ζ ζ ∞ 0 J(ζ, ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ F i (ζ) -F i (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ , (6.1
) whereas the moments of Q ei are expressed as

Q i ei (f ) = i ζ 3 (i -1)f i-2 -(i + 1)f i , i ≥ 1, Q 0 ei = 0. (6.
2)

The expression of Q i ee comes from the linearity of Q ee whereas the moment extraction Q i ei for the electron-ion collision is computed in Appendix C. By setting

f =    f 0 . . . f N    , f =    f 1 . . . f N +1    and Q =    Q 0 . . . Q N    , (6.3)
the moments extraction of 4.3 leads to the system

∂ t f + ζ∂ x f = Q . (6.4)
6.2. Conservation of the realizability domain. We show that the collision operator Q(f ) = Q ee (f ) + Q ei (f ) conserves the realizability domain A. That is why we prove in Proposition 3 that the space homogeneous version of (6.1, 6.3, 6.2, 6.4) preserves A.

Proposition 3. The continuous problem (6.1, 6.3, 6.2, 6.4) considered in the space homogeneous case preserves the realizability domain

A i.e if f (t = 0, ζ) ∈ A then f (t, ζ) ∈ A for any t ≥ 0.
The proof of Proposition 3 is shown by using the following lemma.

Lemma 2. The kinetic equation (4.1, 4.2, 4.3) in the space homogeneous setting preserves the nonnegativity of the solution ie if g(0, ζ, µ) ≥ 0 then g(t, ζ, µ) ≥ 0 for any t ≥ 0.

Proof. For any quantity α, define by α + (resp. α -) its positive part (resp. negative).

Let g be solution to ∂g ∂t = Qee (g) + Q ei (g) , (

with

Qee (g) = 2 3ζ 2 ∂ ζ A(G 0 + , ζ)∂ ζ g + B(G 0 + , ζ)g . (6.6)
Here A and B are given by

A(G 0 + , ζ) = ζ 2 ∞ 0 min( 1 ζ 3 , 1 w 3 )w 4 G 0 (w) + dw , B(G 0 + , ζ) = -ζ 3 ∞ 0 min( 1 ζ 3 , 1 w 3 )w 3 ∂ w g 0 (w) + w 2 dw = 3 ζ 0 w 2 G 0 (w) + dw ,
where G 0 = g 0 /ζ 2 as usual.

For more details on the last relation concerning B we refer ( [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Maxwell-Fokker-Planck equation of electrons[END_REF]). We proceed like in ( [START_REF] Brull | Existence of solutions to discrete coagulation-fragmentation equations with transport and diffusion[END_REF], [START_REF] Mallet | An entropic scheme for an angular moment model for the classical Maxwell-Fokker-Planck equation of electrons[END_REF]). Let H be a convex positive function such as H ′ ≤ 0 and defined by

H(x) = -Cx if x < 0 0 if x ≥ 0 , C > 0 . ( 6
ζ 2 ∂ t H(g)dζ = +∞ 0 ζ 2 Q ee (g)H ′ (g) dζ + +∞ 0 ζ 2 Q ei (g)H ′ (g) dζ . (6.8)
We show firstly the dissipation property

+∞ 0 1 -1 ζ 2 Q ee (g)H ′ (g)dµ dζ ≤ 0. ( 6.9)
The first right-hand side term of (6.8) writes

4π +∞ 0 1 -1 ζ 2 Q ee (g)H ′ (g)dµ dζ = 4π ∞ 0 1 -1 2 3 ∂ ζ A(F 0 + , ζ)∂ ζ g H ′ (g)dµ dζ + +∞ 0 1 -1 2 3 ∂ ζ B(F 0 + , ζ)g H ′ (g)dµ dζ .
(6.10)

Firstly, remark that H can be regularized into a C 2 function satisfying H ′ (0) = H ′′ (0) = 0. In this way, the second derivative of H can be considered.

Hence an integration by parts on the first right-hand side term of (6.10) and the convexity of H lead to

∞ 0 1 -1 2 3 ∂ ζ A(F 0 + , ζ)∂ ζ g H ′ (g)dµ dζ ≤ 0.
Moreover the second term of (6.10) can be simplified as

∞ 0 1 -1 2 3 ∂ ζ B(F 0 + , ζ)g H ′ (g)dµ dζ = ∞ 0 2H ′ (g)∂ ζ g ζ 0 w 2 F 0 + (w)dwdζ + ∞ 0 2H ′ (g)ζ 2 gF 0 + dµ dζ .
An integration by part leads to

+∞ 0 1 -1 2 3 ∂ ζ B(F 0 + , ζ)g H ′ (g)dµ dζ = +∞ 0 1 -1 2ζ 2 F 0 + -H(g) + H ′ (g)g dµ dζ.
Moreover, by construction of H, ∀x, (-H(x) + H ′ (x)x) x = 0. So, (6.9) is satisfied.

Next we show the same inequality as (6.9) for Q ei . By using an integration by parts we get

+∞ 0 1 -1 ζ 2 Q ei (g)H ′ (g) dµ dζ = - +∞ 0 1 -1 ζ 2 (1 -µ 2 ) ∂g ∂µ 2 H ′′ (g)dµdζ .
By convexity of H, Q ei satisfies the inequality (6.9). So

+∞ 0 1 -1 ζ 2 ∂ t H(g) dµ dζ ≤ 0 .
After integrating the previous inequality, we get

+∞ 0 1 -1 ζ 2 H(g(t, ζ)) dµ dζ ≤ +∞ 0 1 -1 ζ 2 H(g(0, ζ)) dµ dζ .
f (0, ζ) = 2πζ 2 1 -1 f (0, ζ, µ)µdµ, f (0, ζ, µ) ≥ 0.
Moreover the solution of (4.1, 4.2, 4.3) for initial condition f (0, ζ, µ) is nonnegative. Next the moment system has a solution belonging to A. Therefore by uniqueness of the solution to (6.4), we get that A is conserved.

6.3. H-theorem. One of the main result of this paper is the following theorem which proves that the system (6.3, 6.4) is entropic.

Theorem 6.1. E = f ln f -f ζ 2 is an entropy for the system (6.3, 6.4). More precisely, we have ∂ t E + ∂ x F ≤ 0 , where F is the entropic flux given by

F = ζ (f ln f -f )µ ζ 2 .
The proof is performed in the same spirit as in ( [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF]). So we refer to this reference for more details.

7. Discretization of the N-moment model. In this section, the moments procedure is applied on the semi-discrete scheme proposed in section 5 when the distribution function is obtained from the minimization entropy principle. So, in the present section we denote for any j ∈ {1; m}, f j by

f j = exp(α j .µ) , α j =    α 0,j . . . α n,j    and α i,j = α i (ζ j ) . (7.1)
7.1. N-moment system for the semi-discretized kinetic equation. By extracting moments on equation (5.1), we obtain the N-moment discretized system

∀j ∈ {1, m} ∂ t f j + ζ j ∂ x f j = Q j , (7.2)
where

f i j = ζ 2 µ i f j f j =    f 0 j . . . f n j    , f j =    f 1 j . . . f n+1 j    , Q j =    Q 0 j . . . Q n j    = µQ j .
Definition 2. The discrete collision operators (Q i ee,j ) j∈{1;m} and (Q i ei,j ) j∈{1;m} involved in (7.2) are given by

       Q i ee,j = G i j+ 1 2 -G i j-1 2 ∆ζj , Q i ei,j = i ζ 3 j ((i -1)f i-2 j -(i + 1)f i j ), (7.3)
where

G i j+ 1 2 = ζ j+ 1 2 m k=0 J(ζ j+ 1 2 , ζ k+ 1 2 ) F 0 k+ 1 2 1 ζ j+ 1 2 F i j+1 -F i j ∆ζ j+ 1 2 - 1 ζ k+ 1 2 F i j+ 1 2 F 0 k+1 -F 0 k ∆ζ k+ 1 2 ζ 2 k+ 1 2
∆ζ k+ 1 2 , (7.4) and

F i j+ 1 2 = 1 -1 f j+ 1 2 µ i dµ . (7.5)
Recall that f j+ 1 2 is defined in (5.4). Remark 5. From a computational point of view, an approximation of (7.5) can be obtained through a usual quadrature formula.

Realizability domain.

In this section, we prove that the discretized collision operator defined in section 5 preserves the realizability domain A. That is why we consider the moment system (7.1, 7.2, 7.3, 7.4, 7.5) in the space homogeneous context. Proposition 4. The space homogeneous version of the semi-discretized problem (7.1, 7.2, 7.3, 7.4, 7.5) preserves the realizability domain A i.e. if f j (t = 0) ∈ A, then f j (t) ∈ A for any t ≥ 0. Lemma 3. The semi-discretized kinetic problem (5.1, 5.2) preserves the nonnegativity of the solution ie if f j (0, µ) ≥ 0 then f j (t, µ) ≥ 0.

Proof. For any j ∈ {1; m}, let g j be solution to

∂g j ∂t = G j+ 1 2 -G j-1 2 ζ 2 j ∆ζ j + Q ei,j , (7.6)
where

G j+ 1 2 = A j+ 1 2 g j -g j+1 ∆ζ j+ 1 2 + B j+ 1 2 g j+1 , B j+ 1 2 = B j+ 1 2 + g + j+ 1 2 g + j+1 , A j+ 1 2 = ζ 2 j+ 1 2 m k=0 min 1 ζ 3 j+ 1 2 , 1 ζ 3 k+ 1 2 F 0 k+ 1 2 + ζ 4 k+ 1 2 ∆ζ k+ 1 2 . (7.7)
Let H be a convex positive function defined as in 6.7. Multiplying equation (7.6) by ζ 2 j H ′ (g j ) summing on j and integrating over µ leads to

1 -1 j ζ j+ 1 2 ∂ t H(g j ) ∆ζ j dµ = 1 -1 j (G j+ 1 2 -G j-1 2 )H ′ (g j )dµ + 1 -1 j ζ 2 j+ 1 2 Q ei,j H ′ (g j ) ∆ζ j+ 1 2 dµ . (7.8)
From (7.7) together to a shift of index, it comes that the first term of the right-hand side of equation (7.8) can be simplified into

1 -1 j (G j+ 1 2 -G j-1 2 )H ′ (g j ) ∆ζ j dµ = 1 -1 j A j+ 1 2 g j -g j+1 ∆ζ j+ 1 2 H ′ (g j+1 ) -H ′ (g j ) dµ + 1 -1 j B j+ 1 2 g j+1 -B j-1 2 g j H ′ (g j ) dµ. (7.9)
The second term of the right hand side of (7.9) can be rewritten as

j B j+ 1 2 g j+1 -B j-1 2 g j H ′ (g j ) = j B j+ 1 2 -B j-1 2 g j H ′ (g j ) + B j+ 1 2 (g j+1 -g j )H ′ (g j ) . (7.10)
Moreover by convexity H, it holds that

H(g j+1 ) ≥ H(g j ) + (g j+1 -g j )H ′ (g j ) .
Then by non-negativity of B j+ 1 2 , equation (7.10) gives

j B j+ 1 2 g j+1 -B j-1 2 g j H ′ (g j ) ≤ j B j+ 1 2 -B j-1 2 g j H ′ (g j ) + B j+ 1 2 (H(g j+1 ) -H(g j ))
. By using a shift of index, we get

j B j+ 1 2 g j+ 1 2 -B j-1 2 g j H ′ (g j ) ≤ j B j+ 1 2 -B j-1 2 (-H(g j ) + g j H ′ (g j ))
.

By construction of H, we get (-H(g j ) + g j H ′ (g j )) = 0 and thus from the previous inequality it comes that

j B j+ 1 2 g j+ 1 2 -B j-1 2 g j H ′ (g j ) ≤ 0 . (7.11)
Moreover since H ′ is a non-decreasing function, the first right hand term side of (7.9) is negative. Therefore, from (7.11), we obtain the following inequality

1 -1 j (G j+ 1 2 -G j-1 2 )H ′ (g j ) ∆ζ j dµ ≤ 0 .
Besides by arguing like for proof of lemma 2, the second term of the right-hand side of equation (7.8) is negative. Then by using definition of H, lemma 3 follows.

Proof. (Proposition 4). We proceed like for the proof of Proposition 3. Consider f j (t = 0) ∈ A. Then

f j (0) = 1 -1 f j (0, µ)µdµ, f j (0, µ) ≥ 0.
From lemma 3, the solution of (5.1, 5.2) for the initial condition f j (0, µ) is nonnegative. So the solution of the moment system belongs to A and we conclude like for the proof of Proposition 3.

7.3. Entropic property of the scheme. Finally, this subsection is devoted to demonstrate the entropic property of the scheme.

Theorem 7.1. E = m j=0 f j ln f j -f j ζ 2 j is an entropy for the system (7.1, 7.

. More precisely, we have

∂ t E + ∂ x F ≤ 0 , where F = m j=0 ζ j (f j ln f j -f j )µ ζ 2 j
is the entropic flux.

Proof. First, after multiplying equation (7.2) by α, we get

α j .∂ t f j + ζ j α j .∂ x f j = α j .Q j .
Hence, by using the expression f j ,

α j .∂ t µf j ζ 2 j + ζ j α j .∂ x µµf j ζ 2 j = α j .µQ j ζ 2 j . (7.12)
The first term of the left hand side term of (7.12) can be expressed as

α j .∂ t µf j = ∂ t α j .µf j -∂ t (α j .µ)f j ,
and simplified into

α j .∂ t µf j = ∂ t log f j f j -∂ t f j . (7.13)
In the same way, the second term of the right hand side can be rewritten into

ζ j α j .∂ x µµf j = ζ j ∂ x ( log f j f j µ ) -ζ j µ∂ x f j . (7.14)
Then, from equations (7.13) and (7.14) and the definition of f j , it comes

∂ t m j=0 f j ln f j -f j ζ 2 j ∆ζ j + ∂ x m j=0 ζ j (f j ln f j -f j )µ ζ 2 j ∆ζ j = log f j Q j ζ 2 j .
From the dissipative property of Q j given in Proposition 2, Theorem 7.1 follows.

7.4. Numerical results. We present in this section a test case, where the use of the M 2 model is relevant compared to the M 1 model. However, using the M 2 model is not so easy than using the M 1 model. Indeed, the computation needs the knowledge of the third order moment f 3 depending of the lower order moments. This relation cannot be computed explicitly and needs some numerical approximations. The Lagrange multipliers (α i ) i=0,2 can be obtained from the moments f i as follows

f (µ) = exp(Σ 2 i=0 α i µ i ) where f i = 1 2 1 µ=-1 µ i f (µ)dµ.
In fact, (α i ) i=0,2 are found by solving the following convex minimization problem:

α = (α i ) i=0,2 = min β 1 µ=-1 exp(Σ 2 i=0 β i µ i )dµ -Σ 2 i=0 β i f i , β = (β i ) i=0,2 . (7.15)
Once the (α i ) i=0,2 have been computed, the moment f 3 can be obtained directly from the undelying distribution function f . Generally, (7.15) is solved by computing the integrals numerically and minimizing the integral by some convex optimization solver (see [START_REF] Hauck | High-Order Entropy-Based Closures for Linear Transport in Slab Geometries[END_REF]). In a first approach, we have used QUADPACK, an adaptative numerical solver from netlib repository (www.netlib.org) and an unconstrained optimization solver HUMSL coming also from netlib. This method will be denoted by M opt N . This method can be very expensive compared to a kinetic computation. Hence we have chosen another approach. Fortunately in the case of M 2 the moments can be computed from the α variables defined in (7.15):

f 0 = 1/2 erf 1/2 -2 α 2 + α 1 √ -α 2 √ πe 1/4 4 α 0 α 2 -α 1 2 α 2 1 √ -α 2 -1/2 erf 1/2 2 α 2 + α 1 √ -α 2 √ πe 1/4 4 α 0 α 2 -α 1 2 α 2 1 √ -α 2 , f 1 = 1/4 e α0 -α1 α 1 √ πe -1/4 α 1 (α 1 -4 α 2 ) α 2 erf 1/2 -2 α 2 + α 1 √ -α 2 (-α 2 ) -3/2 -1/4 e α0 -α1 α 1 √ πe -1/4 α 1 (α 1 -4 α 2 ) α 2 erf 1/2 2 α 2 + α 1 √ -α 2 (-α 2 ) -3/2 -1/2 e α0 -α1 e α2 -e 2 α1 +α2 α 2 , f 2 = 1/8 e α0 -α1 √ πe -1/4 α 1 (α 1 -4 α 2 ) α 2 -2 α 2 + α 1 2 erf 1/2 -2 α 2 + α 1 √ -α 2 (-α 2 ) -5/2 -1/8 e α0 -α1 √ πe -1/4 α 1 (α 1 -4 α 2 ) α 2 -2 α 2 + α 1 2 erf 1/2 2 α 2 + α 1 √ -α 2 (-α 2 ) -5/2
+ 1/4 e α0 -α1 2 e α2 α 2 + α 1 e α2 + 2 e 2 α1 +α2 α 2 -α 1 e 2 α1 +α2 α 2 2 .

In fact f 1 /f 0 , f 2 /f 0 are function of α 1 , α 2 variables on the domain of definition of M 2 . The quantities α 1 , α 2 can be computed on some grids points in the f 1 /f 0 , f 2 /f 0 phase space and then the third reduced moment f 3 /f 0 can be deduced.

Solving the problem at grid points is done by using a two-dimensional dichotomy. The process is rather expensive but used only once. The closure can be obtained by interpolating f 3 /f 0 on the f 1 /f 0 , f 2 /f 0 grid space. This method is called M app 2 . A grid of about 100 points by direction was generated on f 1 /f 0 , f 2 /f 0 phase space.

On the other hand, a same kind of process can be applied for M 3 , leading to f 4 /f 0 interpolation over a 100 3 grid points in the f 1 /f 0 , f 2 /f 0 , f 3 /f 0 phase space. We called this method M app 3 . For the test case, we considered only elastic electon-ion collisions, so we take α ee = 0, α ei = 1. In this case, the model has no ζ derivative. Hence for the sake of simplicity we only consider one point in ζ, ζ = 1. The space interval of study is [0, 5] and we set at x=0, the given incoming boundary conditions

f = 500 exp(-200(1 -µ)), µ ∈ [0, 1].
We have also consider P N approximations which are based on polynomials underlying distribution functions for comparison. These kinds of approximations are very popular in plasma physic. In that case, the P 1 closure is given by f 1 = 1/3f 0 , the P 2 closure by f 3 = 2/5f 1 , and the P 3 closure by f 4 = -3/35f 0 + 6/7f 2 .

For the implementation, we perform a classical Euler scheme in time and we compute the steady state solution of the different equations. The kinetic computation need 180s on a single CPU of i7 Intel Mac computer with GFORTRAN compiler. The moments model M N are rather expensive: 60s for M 1 , 181s for M 2 , and 310s for M 3 . On the other hand, the models M app N are very cheap in comparison (0.7s for M app 1 , 1s for M app 2 , and 2s for M app 3 ) without loss of accuracy. For example, Fig 3 shows the relative error between the two M 2 approaches.

Remark that the most important differences are near the boundaries. To illustrate this fact, the parameters f 1 /f 0 and f 2 /f 0 computed from the kinetic reference solution are plotted on figure 4. We can note that more the solution is far from isotropy (f 1 /f 0 ≈ 0 and f 2 /f 0 ≈ 1/3), more the difference is important. In fact close to the isotropy area, the M 2 closure is a relatively flat function and the interpolation is more accurate than the furthest isotropy area. The distribution function converges to an isotropic distribution function, steady state of the electron-ion collision operator. So the first moment f 1 should be close to 0. However, due to the anisotropic boundary conditions, we observed boundary layer which are represented by both moment models. However, Fig 2 shows that the result given by the M 2 and M 3 model is closer to the kinetic solution than the result given by the M 1 model. 8. Conclusion. Firstly, we have proposed a model for plasma electrons transport with a new consistent collision operator for the electron-electron interactions. Its definition allows mainly to conserve fundamental properties such as the mass and energy conservation and entropy dissipation. Next we have motivated the development of a semi-discretized scheme in the energy variable for the new continuous model. This model still preserves fundamental properties by using a specify entropic average. This M N model is based on an entropy minimization principle and an integration of the kinetic model w.r.t. the angle variable. Finally the discretized moment model constructed from the the semi-discretized kinetic model is proved to be entropic.

Remark that for concrete applications, we prefer low order approach M N models, because they do not use minimization problems which can be too costly.

To extend this model, the forthcoming work should be to consider the model with mobile ions. 9. Acknowledgments. We are thankful with Vladimir Tikhonchuk for fruitful discussions about the pertinence of the physical model.

Hence from (A.17) we get

Φ(V )v = (ζ 2 + ζ ′2 -2ζ ′ v • Ω ′ )v -(ζ 2 v + ζ ′2 (Ω ′ • v)Ω ′ -ζ 2 ζ ′ Ω ′ -ζ ′ (Ω ′ • v)v) (ζ 2 + ζ ′2 -2ζ ′ v • Ω ′ ) 3 2 , = (ζ ′2 -ζ ′ v • Ω ′ )v -(ζ ′2 (Ω ′ • v) -ζ 2 ζ ′ )Ω ′ (ζ 2 + ζ ′2 -2ζ ′ v • Ω ′ ) 3 2 • To achieve the computation of Γ(v) • v, A v = Φ(V )v • v and A Ω ′ = Φ(V )v
• Ω ′ have to be computed. We introduce µ the cosine of that angle between the vector Ω ′ and v, such that Ω ′ • v = ζµ, and we obtain

A v = Φ(V )v • v = (ζ ′2 -ζ ′ v • Ω ′ )ζ 2 + ζ ′ (ζ 2 -ζ ′ (Ω ′ • v))(Ω ′ • v) (ζ 2 + ζ ′2 -2ζ ′ v • Ω ′ ) 3 2 , = ζ ′ (ζ ′ -ζµ)ζ 2 + ζζ ′ µ(ζ 2 -ζζ ′ µ) (ζ 2 + ζ ′2 -2ζζ ′ µ) 3 2 = ζ 2 ζ ′2 (1 -µ 2 ) (ζ 2 + ζ ′2 -2ζζ ′ µ) 3 2
, (A.18)

A Ω ′ = Φ(V )v • Ω ′ = ζ ′ (ζ ′ -v • Ω ′ )(Ω ′ • v) + ζ ′ (ζ 2 -ζ ′ (Ω ′ • v)) (ζ 2 + ζ ′2 -2ζ ′ v • Ω ′ ) 3 2 , = ζ ′ (ζ ′ -ζµ)ζµ + ζ ′ (ζ 2 -ζζ ′ µ) (ζ 2 + ζ ′2 -2ζζ ′ µ) 3 2 = ζ 2 ζ ′ (1 -µ 2 ) (ζ 2 + ζ ′2 -2ζζ ′ µ) 3 2 • (A.19)
Hence by symmetry of Φ(V ) and dΩ ′ = 2πdµ, Γ(v) • v writes

Γ(v) • v = ∞ 0 S 2 Φ(V )v • f (ζ ′ ) v ζ ∂ ζ f (ζ) -f (ζ)Ω ′ ∂ ζ ′ f (ζ ′ ) dΩ ′ ζ ′2 dζ ′ , = ∞ 0 S 2 f (ζ ′ ) 1 ζ ∂ ζ f (ζ)A v -f (ζ)∂ ζ ′ f (ζ ′ )A Ω ′ dΩ ′ ζ ′2 dζ ′ , = ∞ 0 +1 -1 ζ 2 ζ ′2 (1 -µ 2 ) (ζ 2 + ζ ′2 -2ζζ ′ µ) 3 2 f (ζ ′ ) 1 ζ ∂ ζ f (ζ) -f (ζ) 1 ζ ′ ∂ ζ ′ f (ζ ′ ) 2πdµζ ′2 dζ ′ .
As,

+1 -1 ζ 2 ζ ′2 (1 -µ 2 ) (ζ 2 + ζ ′2 -2ζζ ′ µ) 3 2 dµ = 4 3 inf 1 ζ 3 , 1 ζ ′3 , Γ(v) • v can be rewritten as Γ(v) • v = ∞ 0 4π J(ζ, ζ ′ ) f (ζ ′ ) 1 ζ ∂ ζ f (ζ) -f (ζ) 1 ζ ′ ∂ ζ ′ f (ζ ′ ) ζ ′2 dζ ′ . (A.20) As f is isotropic, f 0 (ζ) = ζ 2 S 2 f (ζ) dΩ = 4πζ 2 f (ζ). So ∂ t f 0 = ∂ ζ ζ ∞ 0 J(ζ, ζ ′ ) f 0 (ζ ′ ) ζ ′2 1 ζ ∂ ζ f 0 (ζ) ζ 2 - f 0 (ζ) ζ 2 1 ζ ′ ∂ ζ ′ f 0 (ζ ′ ) ζ ′2 ζ ′2 dζ ′ .
We retrieve exactly the formula for the collision operator involved by the equation on f 0 . We neglect the operator Q ee for the equation on f 1 because we retain only the isotropic part of this operator. The gradient ∇ v f can be rewritten as

∇ v f = ∂f ∂ζ ∇ v ζ + ∂f ∂µ ∇ v µ + ∂f ∂ϕ ∇ v ϕ . (B.21)
Since we have

  ∇ v ζ ∇ v µ ∇ v ϕ   = D(ζ, µ, ϕ) D(v x , v y , v z ) = D(v x , v y , v z ) D(ζ, µ, ϕ) -1 =         µ 1 -µ 2 cos ϕ 1 -µ 2 sin ϕ 1-µ 2 ζ -µ √ 1-µ 2 cos ϕ ζ -µ √ 1-µ 2 sin ϕ ζ 0 -sin ϕ √ 1-µ 2 ζ cos ϕ √ 1-µ 2 ζ        
, the gradient vector can be simplified as Therefore we obtain (4.2) when ∂ ∂ϕ f = 0.

∇ v ζ = -→ Ω , ∇ v µ = 1 -µ 2 ζ -→

R 3

 3

Figure 1 .

 1 Figure 1. Representation of the moments f 0 and f 1 with respect to the energy ζ at the steady state.

1 - 1 ζ 2 H

 112 As g(0, ζ) ≥ 0, the definition of H gives H(g(0, ζ)) = 0. Then by nonnegativity of H, it comes that +∞ 0 (g)dµ dζ = 0. So H(g) = 0 a.e. i.e. g ≥ 0. So G 0+ = G 0 and by uniqueness of the solution of (4.1, 4.2, 4.3) we get the result. Proof. (Proposition 3). Let choose f (0, ζ) ∈ A. Hence

Figure 2 .

 2 Figure 2. Representation of the steady state for the first moment f 1 computed with the M N model, the P N model, for N = 1, 3 and the kinetic one for the space interval [0, 5] and one group in energy, ζ = 1 with 500 points in space and 256 points in µ for the kinetic model.

Figure 3 .

 3 Figure 3. Representation of relative error between M 2 and M app 2

Figure 4 .

 4 Figure 4. Representation of f 1 /f 0 and f 2 /f 0 parameters about the steady state for the kinetic computation.

Appendix B. Expression ( 4 . 2 ) 1 -

 421 of Q ei in spherical coordinate (ζ, µ, ϕ). Consider 2.1 in the homogeneous case with C = C ei . The direction of propagation Ω and the velocity v writeΩ =   µ 1 -µ 2 cos ϕ 1 -µ 2 sin ϕ µ 2 cos ϕ ζ 1 -µ 2 sin ϕ   .

23 ) 3 ∇ 3 Φ 3 ∇ 3 ∇ 3 ∇

 2333333 e µ and ∇ v ϕ = 1 1 -µ 2 ζ -→ e ϕ . (B.22) By using (B.22), (B.21) writes ∇ Next we determine the expression of the eletron/ion operator in its weak form. Consider a test function Ψ. Hence, by using the expression of the divergence term (B.23) and a Green formula, we obtainR v .(Φ(v)∇ v f )Ψ(v)dv = -R v .(Φ(v)∇ v f )Ψ(v)dv =ϕ ). Besides ( -→ e ζ , -→ e µ , -→ e ϕ ) is an orthonormal basis, so (B.24) writes R v .(Φ(v)∇ v f )Ψ(v)dv =µ 2 )ζ 3 dv . (B.25) As dv = ζ 2 dζdµdϕ and by using integration by parts, (B.25) writes R v .(Φ(v)∇ v f )Ψ(v)dv =

  Properties of the collision operator. Next we present the conservation laws and entropy dissipation property of the collision operators C ee and C ei . The electron-electron collision operator satisfies mass, momentum and energy conservation properties

	.4)
	2.2.

supported by NSF grant xx-xxxx.

Appendix A. Equation (3.7). To establish this closure we consider 2.1 in the homogeneous case with the collision operator C = C ee . For the sake of simplicity, the distribution function f is assumed to be isotropic, i.e. f = f (|v|) = f (ζ). Hence we have to deal with

Moreover if the distribution function f is isotropic, Γ(v)v is also isotropic. In order to perform a weak formulation of A.16, consider ϕ some regular test function. Hence, by setting

.

(A.17)

Appendix C. Moment extraction for the electron-ion collision operator.

In this appendix, we show formula 6.2 from (4.2) by using its weak formulation. By considering the i th moment of Q ei and by using successively two Green formula, we get

where ψ denotes some test function.

Appendix D. Proof of Lemma 1. By using the definition (5.6) of B j+ 1 2 , it holds that