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

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

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

An entropic scheme for an angular moment model for the classical Fokker-Planck-Landau equation of electrons

Keywords: entropy minimization, Landau-Fokker-Planck equation, moment systems, entropic scheme

An entropic scheme for an angular moment model for the classical Fokker-Planck equation of electrons

Jessy Mallet, Stéphane Brull, Bruno

Introduction

Classically in kinetic theory, a Fokker-Planck equation is used to describe the evolution of different species in a collisional plasma ( [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF][START_REF] Delcroix | Physique des Plasmas. InterEditions[END_REF]). These charged particles can interact through 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 specie α of particles with velocity v at time t and position x. For the sake of simplicity, we assume in this paper that the plasma consists of electrons and one ion species which are considered fix in the plasma. In order to approximate the solution of such problems, many computational methods have been developed up to now.

There is a variety of models that are used to describe electrons transport: the hydrodynamic ones ( [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]) and the kinetic ones ( [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][START_REF] Duclous | Modélisation et Simulation Numérique multi-échelle du transport cinétique électronique[END_REF][START_REF] Duclous | High order resolution of of the Maxwell-Fokker-Planck-Landau model intended for ICF application[END_REF][START_REF] Kemp | Numerical method for particle simulations at extreme densities and temperatures: weighted particles, relativistic collisions and reduced currents[END_REF]). But in the present paper, we consider an intermediate description beween fluid and kinetic level. Let us consider f 0 the zeroth order moment of the distribution function f (isotropic part), its first order moment f 1 and its second order moment f 2 with respect to the angular part of the velocity variable. The usual model employed in plasma physic can be expressed such as

∂ t f 0 +ζ∂ x f 1 = Q 0 ( f 0 ) , ∂ t f 1 +ζ∂ x f 2 = -κ f 1 , (1.1)
where κ ∈ R + and ζ is the microscopic energy variable. Here Q 0 ( f 0 ) represents an approximation for the electron-electron collision operator. In plasma physics, classical approximation lead to consider that the main contribution for the electron-electron collision operator comes from the isotropic part, that is why here Q 0 depends only of f 0 . This model allows mainly to conserve fundamental properties such as the mass and energy conservation and the entropy dissipation. However this model is ill posed. Indeed whatever the closure chosen, this model does not preserve the realizability domain. Of course, except this model, some kinetic models like Fokker Planck model exist but computations are too much expensive and mainly kinetic codes can not run many collisional times. Most of them are restricted to some decade of collisional time in practice. For example there exists a kinetic code called KETS used in [START_REF] Duclous | Modélisation et Simulation Numérique multi-échelle du transport cinétique électronique[END_REF][START_REF] Duclous | High order resolution of of the Maxwell-Fokker-Planck-Landau model intended for ICF application[END_REF] and cross-validated with an implicit Vlasov-Fokker-Planck code called IMPACT developed by Kingham and Bell [START_REF] Bell | An implicit Vlasov-Fokker-Planck code to model non-local electron transport in 2-D with magnetic fields[END_REF].

Beyond a few decade of collisional time, KETS code becomes too costly. To reduce the computational time, the plasma 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] a bi-fluid compressible Euler model coupled with the Poisson equation is considered. However, for the new high energy target drivers, the kinetic effects are too important to neglect them and replace kinetic calculus by usual macroscopic Euler Models with T i -T e temperatures.

That is why to preserve the realizability domain, we consider in this paper a new approximation of the electron-electron collision operator to correct the model 1.1. This new continuous model satisfies fundamental properties (i.e. conservation laws, entropy dissipation and conservation of the realizability domain). The microscopic velocity is written in spherical coordinates and the model is written by considering moments system with respect to the angular variable. But the choice of the closure on the angle variable is crucial to guarantee reasonable properties for the resulting model. For example, the P 1 model ( [START_REF] Grad | On Kinetic Theory of the Rarefied Gases[END_REF]) neither satisfies positivity of the underlying distribution function of electrons and nor entropy dissipation. Indeed, the closure chosen defines the distribution function by a polynomial function depending on the angle variable. However in [START_REF] Hauck | Positive P N closures[END_REF], a modification of P N model is proposed to correct this defect. Nevertheless the positive P N model does not dissipate the Boltzmann entropy. For this reason there is a considerable motivation to develop an other moment closure. That is why we will derive a new reduced kinetic model called M 1 model whose construction is based on an entropy minimization principle. This model has been firstly introduced in [START_REF] Feugeas | Entropic moment closure hierarchy for the radiative transfert equation[END_REF] to find an approximation of radiative transfer and energy evolution equation. The system of conservation laws that is obtained is closed by minimizing the radiative entropy (see also [START_REF] Frank | Partial moment entropy approximation to radiative transfer[END_REF] for extension). This approach has been extended in [START_REF] Berthon | Numerical methods for balance laws with space dependent flux: application to radiotherapy dose calculation[END_REF] by proposing numerical methods with high order reconstruction providing positivity of the distribution function. In order to get accurate results and costless numerical schemes, a convenient change of variables has been proposed to eliminate the space dependence in the flux function. In kinetic theory, the problem of minimization under moments constraints has been firstly introduced by Levermore in [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]) which could fit transport coefficients at the hydrodynamic limit. Moreover the well-posedness of such problems has been clarified by Junk and Schneider ( [26,[START_REF] Schneider | Entropic approximation in kinetic theory[END_REF]). But in the present case, the aim is different. Here the energy of particles constitutes a free parameter. Then we integrate only the kinetic equation with respect to the angle variable and we return only the energy of particles as kinetic variable. By using a closure defined from a minimization entropy principle, we obtain a closed system.

The main contribution of this paper is firstly the derivation of a new moments system from the Fokker-Planck-Landau system. This model is next shown to be entropic and to preserve the realizability domain. In order to perform numerical simulations, an entropic semi-discretization in the energy variable of the model is proposed. The cornerstone of the construction of the scheme is an entropic average used to define the quantities on the mesh. This definition implies that the Rosenbluth and Landau form of the discretized collision operator are still equivalent. This entropic average has been used for a kinetic equation by Dellacherie in [START_REF] Buet | On the Chang and Cooper scheme applied to a linear Fokker-Planck equation[END_REF][START_REF] Buet | Résolution numérique d'une équation de Fokker-Planck ionique avec température électronique[END_REF][START_REF] Dellacherie | Sur un schéma numerique semi-discret appliqué à un opérateur de Fokker-Planck isotrope[END_REF][START_REF] Dellacherie | Numerical resolution of an ion-electron collision operator in axisymmetrical geometry[END_REF][START_REF] Dellacherie | Numerical solution of an ionic Fokker-Planck equation with electronic temperature[END_REF] when the distribution function is isotropic.

This paper is organized as follows. The principle of the angular moment closure is explained in section 2 on a simple kinetic BGK model. In section 3, the moment model is generalized for the Fokker-Planck-Landau kinetic equation. Here the distribution function is solution of the system constituted by the Fokker-Planck-Landau kinetic equation coupled with the Maxwell equations which describe the time evolution of the electromagnetic field. The system is derived by taking the angular moments defined in section 2. In section 4 the angular moment system is proved to be entropic. Section 5 is devoted to the derivation of an entropic semi-discrete scheme for the space homogeneous moments system. The time is kept as a continuous variable whereas the energy variable is discrete. Next several numerical test cases for the full discrete scheme are performed in section 6. Test cases show that full discrete scheme satisfies the properties of the semidiscrete scheme shown in section 5. Finally we finish in the last section by conclusions and perspectives to this work.

Principle of the angular moment closure

In this section an original derivation of moment closure is proposed. Instead of other existing moment model, the final equation is semi-kinetic. Indeed the system is obtained by an angle variable integration and the system remains kinetic in energy. For the sake of simplicity, we consider a BGK kernel-collision operator to derive the closure. The time evolution of the electron distribution function f (t,x,v), where x ∈ R 3 and v ∈ R 3 , is governed by the following kinetic equation,

∂ t f (t,x,v)+v•∇ x f (t,x,v) = C( f ), (2.1)
where,

C( f ) = ν(|v|)(M f -f ), M f = n m 2πk B T 3/2 exp - mv 2 2k B T , ( 2.2)
and ν(|v|) being the collision frequency.

Next we explain the construction of the M 1 model on the BGK model (2.1, 2.2). This system is constructed by using a minimum entropy principle for angle variable of particle by keeping the energy of particles as kinetic variable. If S 2 is the unit sphere, Ω = v/|v| represents the direction of propagation of the particles. In this case, by setting ζ = |v|, the three first moments with respect to Ω are defined as

f 0 (ζ) = ζ 2 S 2 f (v)dΩ, f 1 (ζ) = ζ 2 S 2 Ω f (v)dΩ, f 2 (ζ) = ζ 2 S 2 Ω⊗Ω f (v)dΩ. (2.3)
Remark that f 0 is a scalar, f 1 is a vector and f 2 is a tensor. By integrating the system (2.1) with respect to Ω it holds that

∂ t f 0 (t,x,ζ) + ∇ x •(ζ f 1 (t,x,ζ)) = ν(ζ)(4πζ 2 M f (t,x,ζ)-f 0 (t,x,ζ)), ∂ t f 1 (t,x,ζ) + ∇ x •(ζ f 2 (t,x,ζ)) = -ν(ζ) f 1 (t,x,ζ).
Of course, the open problem remains the computation of moment f 2 as a function of f 0 and f 1 . This problem is overcome by solving the entropy minimisation problem

min g≥0 H(g)/∀ζ ∈ R + , S 2 g(Ω,ζ)dΩ = f 0 (ζ), S 2 Ω g(Ω,ζ)dΩ = f 1 (ζ) , (2.4)
where

H( f ) = R 3 ( f ln( f )-f ) dv (2.5)
represents the Boltzmann entropy. Hence the minimum entropy principle [START_REF] Levermore | Moment Closure Hierarchies for Kinetic Theories[END_REF][START_REF] Minerbo | Maximum entropy Eddington Factors[END_REF] applied to (2.4) implies that f should be taken under the following form,

f = ρ(ζ)exp(-Ω•α(ζ)), ( 2.6)
where ρ is a non negative scalar (ρ≥0), and α is a three component real valued vector. An important parameter is the anisotropy parameter a = f 1 / f 0 , which satisfies by construction |a| ≤ 1. The computation of the different moments of the distribution function (2.6) gives

f 0 = 4πρ sinh(|α|) |α| , f 1 = 4πρ sinh(|α|)(1-|α|coth(|α|)) |α| 3 α.
(2.7)

Figure 1: χ given by formula (2.10).

In fact, these relations can be combined to give,

a = 1-|α|coth(|α|) |α| 2 α, ( 2.8)
or by taking the modulus of a,

|a| = |α|coth(|α|)-1 |α| . (2.9)
The relation (2.9) cannot be inverted explicitly by hand. However, this relation determines a unique solution which can be computed and tabulated numerically. a is assumed to be known, f 2 can also be computed

f 2 = f 0 1-χ 2 I + 3χ-1 2
f 1 | f 1 | ⊗ f 1 | f 1 | , where χ = |α| 2 -2|α|coth(|α|)+2 |α| 2 .
(2.10) χ factor is plotted on figure 1.

Remark 2.1. Notice that if χ = 1/3, we get the P 1 model ( [START_REF] Grad | On Kinetic Theory of the Rarefied Gases[END_REF]) described in introduction.

Moment model for the classical case

In this section we recall some backgrounds about the classical kinetic models for plasma physic. Firstly, the Maxwell-Fokker-Planck-Landau model is presented. In particular, the stationary states of the collision are precised. In a second subsection the M 1 model associated to Fokker-Planck-Landau is presented.

Classical Fokker-Planck-Landau kinetic equation

The time evolution of the electron distribution function f (t,x,v) is governed by the Fokker-Planck-Landau kinetic equation which writes ( [8, 14, 28])

∂ t f +v•∇ x f = C ee ( f , f )+C ei ( f ), (3.1)
where C ee and C ei are the electron-electron and electron-ion collision operators. The equation (3.1) can be coupled to the Maxwell equations through the Lorentz force. The electron-electron collision operator is given by

C ee ( f , f ) = α ee ∇ v • v ′ ∈R 3 |u| 2 I -u⊗u |u| 3 f (v ′ )∇ v f (v)-f (v)∇ v ′ f (v ′ ) dv ′ , (3.2)
where u = v ′ -v is the relative velocity of electrons, I is the unit tensor and the coefficient α ee is a nonnegative physical constant.

The second operator is related to the electron-ion interactions,

C ei ( f ) = α ei ∇ v • |v| 2 I -v⊗v |v| 3 ∇ v f (v) , (3.3)
where the coefficient α ei is also a nonegative.

Next we present the fundamental properties of the collision operator C ee and C ei (see [START_REF] Buet | Numerical analysis of conservative and entropy schemes for the Fokker-Planck-Landau equation[END_REF][START_REF] Degond | An entropy scheme for the fokker-planck collision operator of plasma kinetic theory[END_REF] for more details).

The electron-electron collision operator satisfies mass, momentum and energy conservations property

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

R 3 C ei ( f ) 1 |v| 2 dv = 0, t ≥ 0 .
They both satisfy the entropy dissipation property, i.e.

R 3 C ei ( f ) log f dv ≤ 0 and R 3 C ee ( f , f ) log f dv ≤ 0 .
Therefore H( f ) defined in (2.5) is a Lyapounov function for (3.1, 3.2, 3.3).

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

f = n m e 2πk B T 3 2 exp - m e (v-u e ) 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. The equilibrium state of the electron-ion collision operator C ei (i.e. C ei ( f ) = 0) is given by the isotropic functions f = f (|v|) .

The equilibrium state considered by both collision operators is given by the isotropic Maxwellian distribution function

f = n m e 2πk B T 3 2 exp - m e v 2 2k B T .

Moment closure for the classical kinetic equation

As the main part of angular closure is given in section 2, we give only the final result of the closure on Fokker-Planck-Landau system (3.1, 3.2, 3.3). For f 0 , f 1 and f 2 defined as in (2.3), we obtain

∂ t f 0 +∇ x •(ζ f 1 ) = Q 0 ( f 0 ), (3.4) ∂ t f 1 +∇ x •(ζ f 2 ) = Q 0 ( f 1 )+Q 1 ( f 1 ), (3.5)
where the collision operators Q 0 and Q 1 are given by,

Q 0 ( f i ) = 2α ee 3 ∂ ζ ζ 2 A(ζ)∂ ζ f i ζ 2 -ζB(ζ) f i , i = 0,1, (3.6
)

Q 1 ( f 1 ) = - 2α ei ζ 3 f 1 . (3.7)
The coefficients α ee and α ei are introduced in the formula (??) and (??), and the coefficients A(ζ) and B(ζ) are given by

A(ζ) = ∞ 0 min( 1 ζ 3 , 1 w 3 )w 2 f 0 (w)dw, (3.8)
B(ζ) = ∞ 0 min( 1 ζ 3 , 1 w 3 )w 3 ∂ w f 0 (w) w 2 dw. (3.9)
The expression of coefficient B(ζ) (see appendix A) can be simplified as

B(ζ) = - 3 ζ 3 ζ 0 f 0 (w)dw.
The equation (3.4) represents the balance of the charge and involves the electron-electron collision kernel C ee , whereas the equation (3.5) represents the balance of the momentum and involves the electron-ion collision kernel C ei .

Entropic properties of the continuous model

Firstly we prove that the realizability domain is preserved by the continuous model (3.4, 3.5, 3.6, 3.7). Next the model is shown to satisfy a H theorem.

Introduction

From now, we chose a one dimensional direction of propagation. Therefore we take µ ∈ [-1,1] as the direction of propagation variable instead of Ω. In that case the formula (2.6) writes

f = ρ exp(-αµ), ρ ∈ R + , α ∈ R, µ ∈ [-1,1], (4.1)
and the closure relation writes f 2 = f 0 χ. Next, we set

F i (ζ) = f i (ζ) ζ 2 = 1 -1 µ i ρ(ζ)exp(-α(ζ)µ)dµ. (4.2)
By using the expressions of (3.8, 3.9) into (3.6) leads to the following new expressions for

Q 0 ( f i ), i ∈ {1,2}.            Q 0 ( f 0 ) = ∂ ζ ζ ∞ 0 J(ζ,ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ (F 0 (ζ))-F 0 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ , Q 0 ( f 1 ) = ∂ ζ ζ ∞ 0 J(ζ,ζ ′ ) F 0 (ζ ′ ) 1 ζ ∂ ζ (F 1 (ζ))-F 1 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ζ ′2 dζ ′ , (4.3)
with

J(ζ,ζ ′ ) = 2α ee 3 inf 1 ζ 3 , 1 ζ ′3 ζ ′2 ζ 2 . (4.4)

Realizability domain

Firstly, introduce the realizability domain A defined by

A = F = F 0 F 1 ∈ R 2 /∃ f : [-1,1] → R + such that F i = µ i f (α,µ) ∀i = 0,1 .
Remark 4.1. In the particular case of moments up to order 1, the realizability domain is shown to be equal to B defined by

B = F = F 0 F 1 ∈ R 2 , F 0 > 0 and |F 1 | < F 0 ∪{(0,0)} . Proposition 4.1. the realizability domain A is conserved by the system (3.4, 3.5, 3.6, 3.7) i.e. if ( f 0 (0,ζ), f 1 (0,ζ)) ∈ A, then ( f 0 (t,ζ), f 1 (t,ζ)) ∈ A, ∀t .
In order to achieve the proof of the proposition 4.1, we aim to prove the following result.

Lemma 4.1. Let g be a solution to the problem

∂g ∂t = Q(g), (4.5) with Q(g) = 2α ee 3 ∂ ζ ζ 2 A( f 0 ,ζ)∂ ζ g ζ 2 -ζB( f 0 ,ζ)g , (4.6)
where A( f 0 ,ζ) and B( f 0 ,ζ) are defined in (3.8, 3.9). Therefore if g(0) ≥ 0, then g(t) ≥ 0 ∀t.

Proof of Lemma. 4.1. By setting G = g ζ 2 and A(F 0 ,ζ) = ∞ 0 min( 1 ζ 3 , 1 w 3 )w 4 F 0 (w)dw and B(F 0 ,ζ) = 3 ζ 0 w 2 F 0 (w)dw, equations (4.5, 4.6) leads to ζ 2 ∂ t G = 2α ee 3 ∂ ζ ζ 2 A(F 0 ,ζ)∂ ζ G+ B(F 0 ,ζ)G . (4.7)
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 ∈ R + .
By multiplying equation (4.7) by H ′ (G) and integrating over ζ, we get

∞ 0 ζ 2 ∂ t H(G)dζ = ∞ 0 2α ee 3 ∂ ζ ζ 2 A(F 0 ,ζ)∂ ζ G H ′ (G)dζ + ∞ 0 2α ee 3 ∂ ζ B(F 0 ,ζ)G H ′ (G)dζ. (4.8)
By using integration by parts on the first right-hand side term of (4.8) together with the convexity of H implies the inequality

∞ 0 2α ee 3 ∂ ζ ζ 2 A(F 0 ,ζ)∂ ζ G(ζ) H ′ (G)dζ ≤ 0.
Moreover the second term of (4.8) can be simplified into

∞ 0 2α ee 3 ∂ ζ B(F 0 ,ζ)G H ′ (G)dζ = ∞ 0 2α ee H ′ (G)∂ ζ G ζ 0 w 2 F 0 (w)dwdζ + ∞ 0 2α ee ζ 2 H ′ (g)ζ 2 GF 0 dζ .
By using an integration by part, we get

∞ 0 2α ee 3 ∂ ζ B(F 0 ,ζ)G H ′ (G)dζ = ∞ 0 2α ee ζ 2 G -H(G)+ H ′ (G)G dζ.
Moreover, by construction of H, ∀x, (-

H(x)+ H ′ (x)x)x = 0. Therefore, (4.8) gives ∞ 0 ζ 2 ∂ t H(G)dζ ≤ 0 .
After integrating the previous inequality, we get

∞ 0 ζ 2 H(G(t))dζ ≤ ∞ 0 ζ 2 H(G(t = 0))dζ .
By using the definition of H,

∞ 0 H(G)dζ = 0.
Then H(G) ≥ 0 a.e. and the positivity of G is preserved.

Let us now finish the proof of proposition 4.1:

Proof. By using lemma 4.1 with g = f 0 , the positivity of F 0 follows. For the sake of clarity, let us denote g = f 0 + f 1 . By adding equations (3.4, 3.5), it comes that

∂g ∂t = 2α ee 3 ∂ ζ ζ 2 A( f 0 ,ζ)∂ ζ g ζ 2 -ζB( f 0 ,ζ)g - 2α ei ζ 3 f 1 . (4.9)
By proceeding as in the proof of lemma 4.1, we get

∞ 0 ζ 2 ∂ t H(G)dζ ≤ - ∞ 0 2α ei ζ F 1 H ′ (G)dζ . Moreover, if F 0 +F 1 ≥ 0 then H ′ (G) = 0. F 0 +F 1 ≤ 0, hence we get F 1 ≤ 0 and as H ′ ≤ 0, we obtain finally ∞ 0 ζ 2 ∂ t H(G)dζ ≤ 0 .
Therefore by arguing as previously, we get that the positivity of G is preserved. By proceeding in the same way with g= f 0f 1 and g= f 0 + f 1 , it comes that |F 1 |≤F 0 . But to end proof of the conservation of the realizability domain, it remains to show that the last inequality is a strict inequality. Hence, we consider g = f 0 + f 1 such that there exists (t 0 ,ζ 0 ) such that g(t 0 ,ζ 0 ) = 0 . We have shown previously that g ≥ 0. So we get

∂g ∂ζ (t 0 ,ζ 0 ) = 0, ∂ 2 g ∂ζ 2 (t 0 ,ζ 0 ) ≥ 0 , ∂g ∂t (t 0 ,ζ 0 ) = 0 .
Therefore equation (4.9) can be simplified into

2α ee 3 ζ 2 A( f 0 ,ζ)∂ 2 ζ g ζ 2 + 2α ei ζ 3 f 0 = 0 . (4.10)
Since each term of the left-hand side of equation (4.10) is positive, it comes that f 0 =0. So f 1 =0. This situation corresponds to (0,0) which belongs to the realizability domain. The case g = f 0f 1 can be analogously treated.

H Theorem

The main result of this subsection is the following Theorem proving that (3.4, 3.5, 4.3) is entropic.

Theorem 4.1. E = ln(ρ) f 0 -α f 1 is an entropy for the system (3.4, 3.5, 4.3). More precisely, we have

∂ t E+∇•F ≤ 0,
where F is the entropic flux given by F = ln(ρ)

f 1 -α f 2 .
Proof. Multiply equation (3.4) by ln(ρ) and equation (3.5) by -α and add the two equalities. In that case, the left-hand side is equal to ∂ t E+∇•F and it remains to prove that

+∞ 0 ln(ρ)Q 0 ( f 0 )dζ - +∞ 0 α(ζ)Q 0 ( f 1 )dζ - +∞ 0 α(ζ)Q 1 ( f 1 )dζ ≤ 0. (4.11)
By using the definition of Q 1 ( f 1 ), equation (2.8) and that α f 1 ≤ 0, it comes that

- +∞ 0 α(ζ)Q 1 ( f 1 ))dζ ≤ 0.
Next to prove that

+∞ 0 ln(ρ)Q 0 ( f 0 )dζ ≤ 0, we use a Green formula in the expression of +∞ 0 ln(ρ)Q 0 ( f 0 )dζ and we get +∞ 0 ∂ ζ ζ +∞ 0 J(ζ,ζ ′ ) f 0 (ζ ′ ) ζ ′2 1 ζ ∂ ζ ( f 0 (ζ) ζ 2 )- f 0 (ζ) ζ 2 1 ζ ′ ∂ ζ ′ ( f 0 (ζ ′ ) ζ ′2 ) (ζ ′ ) 2 dζ ′ ln(ρ)dζ = - +∞ 0 +∞ 0 J(ζ,ζ ′ ) 1 ζ F 0 (ζ ′ )∂ ζ F 0 (ζ)- 1 ζ ′ F 0 (ζ)∂ ζ ′ F 0 (ζ ′ ) ∂ ζ ρ ρ 1 ζ (ζ ′ ) 2 dζdζ ′ . (4.12) Next we compute 1 ζ F 0 (ζ ′ )∂ ζ F 0 (ζ)- 1 ζ F 0 (ζ ′ )∂ ζ F 0 (ζ). From (4.
2), we get the relation

∂ ζ F 0 (ζ) = 1 -1 ∂ ζ ρ(ζ)exp(-α(ζ)µ)dµ- 1 -1 µ∂ ζ α(ζ)ρ(ζ)exp(-α(ζ)µ)dµ. (4.13)
The expressions of F 0 and

∂ ζ F 0 give 1 ζ F 0 (ζ ′ )∂ ζ F 0 (ζ)- 1 ζ F 0 (ζ)∂ ζ ′ F 0 (ζ ′ ) = 1 -1 1 -1 ρ(ζ)exp(-α(ζ)µ)ρ(ζ ′ )exp(-α(ζ ′ )µ ′ ) ∂ ζ ρ(ζ) ρ(ζ)ζ - µ ζ ∂ ζ α(ζ)- ∂ ζ ′ ρ(ζ ′ ) ρ(ζ ′ )ζ ′ + µ ′ ζ ′ ∂ ζ ′ α(ζ ′ ) dµdµ ′ . (4.14)
Next by setting

K(ζ,ζ ′ ,µ,µ ′ ) = J(ζ,ζ ′ )ζ 2 ζ ′2 ρ(ζ)exp(-α(ζ)µ)ρ(ζ ′ )exp(-α(ζ ′ )µ ′ ), (4.15) δ(ζ) = ∂ ζ ρ(ζ) ρ(ζ) 1 ζ , β(ζ) = ∂ ζ α(ζ) ζ . (4.16)
and by using equality (4.14) in (4.12) we get

- +∞ 0 ζ +∞ 0 J(ζ,ζ ′ )( 1 ζ F 0 (ζ ′ )∂ ζ F 0 (ζ)- 1 ζ ′ F 0 (ζ)∂ ζ ′ F 0 (ζ ′ )) ∂ ζ ρ ρζ dζdζ ′ = - +∞ 0 +∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) δ(ζ)-δ(ζ ′ δ(ζ)dζdζ ′ dµdµ ′ + +∞ 0 +∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) µβ(-ζ)-µ ′ β(ζ ′ ) δ(ζ)dζdζ ′ dµdµ ′ . The transform (ζ,ζ ′ ) → (ζ ′ ,ζ) leads to - +∞ 0 ζ +∞ 0 J(ζ,ζ ′ )( 1 ζ F 0 (ζ ′ )∂ ζ F 0 (ζ)- 1 ζ ′ F 0 (ζ)∂ ζ ′ F 0 (ζ ′ )) ∂ ζ ρ ρζ dζdζ ′ = - 1 2 +∞ 0 +∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) δ(ζ)-δ(ζ ′ ) 2 dζdζ ′ dµdµ ′ + 1 2 +∞ 0 +∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) µβ(ζ)-µ ′ β(ζ ′ ) δ(ζ)-δ(ζ ′ ) dζdζ ′ dµdµ ′ .
(4.17)

Next, for the remaining term

+∞ 0 Q 0 ( f 1 )(-α(ζ))dζ = +∞ 0 +∞ 0 J(ζ,ζ ′ )( 1 ζ F 0 (ζ ′ )∂ ζ F 1 (ζ)- 1 ζ ′ F 1 (ζ)∂ ζ ′ F 0 (ζ ′ )) ∂ ζ α ζ (ζ ′ ) 2 dζdζ ′ ,
we proceed as previously. The expression of F 1 given in (4.2) leads to

∂ ζ F 1 (ζ) = 1 -1 µ∂ ζ ρ(ζ)exp(-α(ζ)µ)dµ- 1 -1 µ 2 ∂ ζ α(ζ)ρ(ζ)exp(-α(ζ)µ)dµ. (4.18)
Therefore by using expressions (4.13, 4.18), we get

+∞ 0 ∞ 0 J(ζ,ζ ′ )ζ 2 ζ ′2 F 0 (ζ ′ ) 1 ζ ∂ ζ F 1 (ζ)-F 1 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ∂ ζ α ζ dζdζ ′ = +∞ 0 ∞ 0 1 -1 1 -1 2α ee 3 K(ζ,ζ ′ ,µ,µ ′ ) δ(ζ)-δ(ζ ′ ) µβ(ζ)dζdζ ′ dµdµ ′ + +∞ 0 ∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) β(ζ ′ )µ ′ -β(ζ)µ µβ(ζ)dζdζ ′ dµdµ ′ .
Then by using again the transform

(ζ,ζ ′ ) → (ζ ′ ,ζ) it comes that +∞ 0 ∞ 0 J(ζ,ζ ′ )ζ 2 ζ ′2 F 0 (ζ ′ ) 1 ζ ∂ ζ F 1 (ζ)-F 1 (ζ) 1 ζ ′ ∂ ζ ′ F 0 (ζ ′ ) ∂ ζ α ζ dζdζ ′ = 1 2 +∞ 0 ∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) δ(ζ)-δ(ζ ′ ) µβ(ζ)-µ ′ β(ζ ′ ) dζdζ ′ dµdµ ′ - 1 2 +∞ 0 ∞ 0 1 -1 1 -1 K(ζ,ζ ′ ,µ,µ ′ ) β(ζ ′ )µ ′ -β(ζ)µ 2 dζdζ ′ dµdµ ′ . (4.19)
Finally, we add the right-hand sides of (4.17) and (4.19) and by using the inequality

(δ(ζ)-δ(ζ ′ ))(β(ζ)µ-β(ζ ′ )µ ′ ) ≤ 1 2 ((δ(ζ)-δ(ζ ′ )) 2 +(β(ζ)µ-β(ζ ′ )µ ′ ) 2 ), (4.20)
we obtain the non positivity of (4.11) and the theorem follows.

Semi-discretization of the problem

In this section we propose an energy discretization such that the Landau form and the Rozenbluth form of the discretized operators are the same (see Remark 5.1). The key point of the scheme construction is the approximation of F 0 and F 1 on the dual mesh through an entropic average. .

Introduction

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

1 2 =]ζ i-1 ,ζ i [, ∀i ∈ [1,m]
, where ζ i = i∆ζ i and m ∈ N the number of points which discretize the energy domain. ∆ζ i represents the step of discretization which can be variable. We denote by D its associated dual mesh consisting of cells ). For the sake of clarity, we define

D i =]ζ i-1 2 ,ζ i+ 1 2 [ where ζ i-1 2 = (i-
F k i (ζ) = f k i (ζ) ζ 2 i , k = 0,1. (5.1)
The aim of this section is to propose an entropic discretization of (3.4, 3.5, 3.6, 3.7)

   ∂ t f 0 i = Q 0 i , ∂ t f 1 i = Q 1 i , (5.2)
where the operators (4.3). The expressions of Q 0 i and Q 1 i are given in (5.3, 5.4).

Q 0 i and Q 1 i are respectively the discrete forms of Q 0 ( f 0 ) and Q 0 ( f 1 )+ Q 1 ( f 1 ) defined in
Definition 5.1. The collision operators Q 0 i and

Q 1 i involved in (5.2) are respectively given by                Q 0 i = G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i , Q 1 i = G 1 i+ 1 2 -G 1 i-1 2 ∆ζ i - 2α ei ζ 3 i f 1 i , (5.3) with                    G 0 i+ 1 2 = ζ i+ 1 2 ∑ j J(ζ i+ 1 2 ,ζ ′ j+ 1 2 )   F 0 j+ 1 2 1 ζ i+ 1 2 F 0 i+1 -F 0 i ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 F 0 i+ 1 2 F 0 j+1 -F 0 j ∆ζ ′ j+ 1 2   ζ ′2 j+ 1 2 ∆ζ ′ j+ 1 2 , G 1 i+ 1 2 = ζ i+ 1 2 ∑ j J(ζ i+ 1 2 ,ζ ′ j+ 1 2 )   F 0 j+ 1 2 1 ζ i+ 1 2 F 1 i+1 -F 1 i ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 F 1 i+ 1 2 F 0 j+1 -F 0 j ∆ζ ′ j+ 1 2   ζ ′2 j+ 1 2 ∆ζ ′ j+ 1 2 , (5.4)
where

F k j+ 1 2 = µ k f j+ 1 2 (µ) , k ∈ {0;1}, f j+ 1 2 (µ) = f j+1 (µ)-f j (µ) ln( f j+1 )(µ)-ln( f j )(µ) , (5.5
)

J(ζ i+ 1 2 ,ζ ′ j+ 1 2 ) = 2α ee 3 inf   1 ζ 3 i+ 1 2 , 1 ζ ′3 j+ 1 2   ζ 2 i+ 1 2 ζ ′2 j+ 1 2 .
Remark 5.1. One important point is the following. When f is isotropic, the collision operator Q 0 ( f ) can be written with its Landau form

∂ ζ ζ ′2 ∞ 0 J(ζ,ζ ′ ) f (ζ ′ ) ζ ′ ∂ ζ f (ζ) ζ 2 - f (ζ) ζ ∂ ζ ′ f (ζ ′ ) ζ ′2 ,
or in its Rosenbluth form as

∂ ζ ζ ′2 ∞ 0 J(ζ,ζ ′ ) f (ζ) f (ζ ′ ) ∂ ζ ln f (ζ) ζ 2 -∂ ζ ′ ln f (ζ ′ ) ζ ′2 .
These relations are of course equivalent for the continuous model, but not always for the discrete one. This problem has already been first considered in ( [START_REF] Buet | Résolution numérique d'une équation de Fokker-Planck ionique avec température électronique[END_REF][START_REF] Dellacherie | Numerical solution of an ionic Fokker-Planck equation with electronic temperature[END_REF]) by using an entropic average (see extensions in [START_REF] Buet | On the Chang and Cooper scheme applied to a linear Fokker-Planck equation[END_REF][START_REF] Dellacherie | Sur un schéma numerique semi-discret appliqué à un opérateur de Fokker-Planck isotrope[END_REF][START_REF] Dellacherie | Numerical resolution of an ion-electron collision operator in axisymmetrical geometry[END_REF]). In the present paper, the entropic average is taken on the underlying distribution function.

Realizability domain

Proposition 5.1. The realizability domain is conserved with respect to time by the semidiscretized problem (5.1, 5.2, 5.3, 5.4).

Proof. Let F 0 such as F 0 i (0) >0 for any i. Hence there is τ >0, s.t. F 0 i (t) >0 for any t ∈ [0,τ[ and any i. Suppose that there exist i 0 such that F 0 i 0 (τ)=0 and assume that i 0 is the smallest index i s.t. F 0 i (τ) = 0. By computing the time evolution of f 0 i in ζ i 0 , we get

∂ t f 0 i 0 (τ) = 1 ∆ζ A i 0 + 1 2 F 0 i 0 +1 (τ)-F 0 i 0 (τ) ∆ζ -A i 0 -1 2 F 0 i 0 (τ)-F 0 i 0 -1 (τ) ∆ζ + 1 ∆ζ B i 0 + 1 2 F 0 i 0 + 1 2 (τ)-B i 0 -1 2 F 0 i 0 -1 2 (τ) , (5.6)
where

A i 0 + 1 2 = ∑ j J(ζ i 0 + 1 2 ,ζ ′ j+ 1 2 )F 0 j+ 1 2 (τ)∆ζ ′ j+ 1 2 ,
and

B i 0 + 1 2 = ζ i 0 + 1 2 ∑ j J(ζ i 0 + 1 2 ,ζ ′ j+ 1 2 ) F 0 j+1 (τ)-F 0 j (τ) ∆ζ ′ j+ 1 2 ∆ζ ′ j+ 1 2 .
From the definition of

F 0 i 0 , if F 0 i 0 → 0, then f i 0 → 0. So from (5.5), it holds that f i 0 + 1 2 ≤ f i 0 . Hence F i 0 + 1 2
→ 0. Therefore, (5.6) can be rewritten as

∂ t f i 0 (τ) = 1 ∆ζ 2 A i 0 + 1 2 F 0 i 0 +1 (τ)+ A i 0 -1 2 F 0 i 0 -1 (τ) . ( 5

.7)

As A i 0 + 1 2 ≥0 and A i 0 -1 2 ≥0, equation (5.7) gives ∂ t f i 0 (τ)≥0. F 0 satisfy for ∆t small enough

F 0 i 0 (τ-∆t) = F 0 i 0 (τ)-∆t∂ t F 0 i 0 (τ)+O(∆t 2 ) = -∆t∂ t F 0 i 0 (τ)+O(∆t 2 ) , because F 0 i 0 (τ) = 0.
If we consider that ∂ t f i 0 (τ) > 0, we get F 0 i 0 (τ-∆t) < 0 for small ∆t. Therefore we have

1 ∆ζ 2 A i 0 + 1 2 F 0 i 0 +1 (τ)+ A i 0 -1 2 F 0 i 0 -1 (τ) = 0. If A i 0 + 1 2 = 0 or A i 0 -1 2 = 0, the definition of A i 0 + 1 2
gives ∀i F 0 i (τ) = 0. In the other case, we get F i 0 +1 (τ) = F i 0 -1 (τ) = 0. By induction on i we still obtain that ∀i F i (τ) = 0. Therefore by using the Cauchy-Lipschitz theorem to the Cauchy problem (5.7) with the initial conditions ∀i, F i (τ) = 0, we get that F i (t) = 0, ∀i,∀t > τ.

If for any i F i (0) = 0, then F i (t) = 0 is a global solution. Finally consider the case where F i 0 (0)=0 and F i 0 +1 (0)>0. If there is τ>0 s.t. F i 0 ≥0 on ]0,τ], we are in the same situation as the first one. If F i 0 <0 on ]0,τ], we have for ∆t small enough, F i 0 (∆t)=∂ t F i 0 (0)∆t+O(∆t 2 ). Then ∂ t F i 0 (0) ≤ 0. As ∂ t F i 0 (0) < 0, is impossible because of equation (5.7), ∂ t F i 0 (0) = 0. But this situation leads to F i 0 +1 (0) =0 which contradicts our assumption. The case where ∂ t F i 0 -1 (0) < 0 and ∂ t F i 0 (0) = 0 being treated analogously, the positivity of the initial condition is preserved for f 0 . By arguing like for the continuous model we conclude that the realizability domain is preserved.

Boundary conditions

We have shown in the previous subsection that the semi-discretized model (5.1, 5.2, 5.3, 5.4) conserves the realizability domain. So we have to choose boundary conditions adapted to the conservation of the realizability domain. Consider here the space homogeneous system. The conservation laws for electron-electron collision operator (3.4) read in spherical coordinates

S 2 ∞ 0 C ee ( f , f )     1 ζΩ ζ 2 2     dΩζ 2 dζ = 0, t ≥ 0 .
(5.8)

Therefore, the conservations of mass, momentum and energy write

∞ 0 Q 0 ( f 0 ) 1 ζ 2 2 ζ 2 dζ = 0 and ∞ 0 Q 0 ( f 1 )ζζ 2 dζ = 0,
which can be rewritten in discretized form as

∞ ∑ j=0 Q 0 ( f 0 j ) 1 ζ 2 j 2 ζ 2 j ∆ζ j = 0 and ∞ ∑ j=0 Q 0 ( f 1 j )ζ j ζ 2 j ∆ζ j = 0.
After integrating over ζ j the first equation of (5.2), mass conservation imposes the boundary condition G 0

-1 2 = 0.
After multiplying by ζ j and integrating over ζ j the second equation of (5.2), we obtain

∞ ∑ j=0 ζ 3 j ∂F 1 j ∂t ∆ζ j = ∞ ∑ j=0 ζ j Q 1 j ∆ζ j = G 1 -1 2 .
However, to satisfy the proof of Proposition 4.1 on the realizability domain, we must have same conditions on G 0 and G 1 . So we use an approximation for f 1 , i.e. we take

G 1 -1 2 =0.
Of course, we lost momentum conservation but this conservation is not satisfied by the operator C ei and so by C ee +C ei .

Entropic property of the scheme

Theorem 5.1. For G defined in formula (5.4),

E = ∑ i f i ln( f i )-f i ζ 2 i = ∑ i G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i ∆ζ i ln(ρ i )-∑ i G 1 i+ 1 2 -G 1 i-1 2 ∆ζ i ∆ζ i α i + 2α ei ζ 3 i α i f 1 i (5.9)
is an entropy for the system (5.2, 5.3, 5.4).

Proof. Firstly the third term defining E in (5.9) satisfies

2α ei ζ 3 i α i f 1 i ≤ 0 .
Next the first term in (5.9) can be rewritten

∑ i G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i ∆ζ i ln(ρ i ) = -∑ i G 0 i+ 1 2 ln ρ i+1 ρ i .
(5.10)

By using the discrete form of the collision operator electron-electron G 0 ( f 0 ) given by (5.4), the right-hand side of (5.10) reads

-∑ i,j ζ i+ 1 2 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 )   F 0 j+ 1 2 1 ζ i+ 1 2 F 0 i+1 -F 0 i ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 F 0 i+ 1 2 F 0 j+1 -F 0 j ∆ζ ′ j+ 1 2   ln ρ i+1 ρ i ζ ′2 j+ 1 2 ∆ζ ′ j+ 1 2 . (5.11)
By using the expression of F i+ 1 2

, we obtain

∑ i G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i ∆ζ i ln(ρ i ) = -∑ i,j 1 -1 1 -1 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 )ζ 2 i+ 1 2 ζ ′2 j+ 1 2 f i+ 1 2 (µ) f j+ 1 2 (µ ′ ) 1 ζ i+ 1 2 f i+1 (µ)-f i (µ) f i+ 1 2 (µ ′ )∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 f j+1 (µ)-f j (µ ′ ) f j+ 1 2 (µ ′ )∆ζ ′ i+ 1 2 1 ζ i+ 1 2 (ln(ρ i+1 )-ln(ρ i )) ∆ζ i+ 1 2 ∆ζ i+ 1 2 ∆ζ ′ j+ 1 2 dµdµ ′ .
According to the definition of f j+ 1 2

given in (5.5), we get

∑ i G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i ∆ζ i ln(ρ i ) = -∑ i,j 1 -1 1 -1 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 ) f i+ 1 2 (µ) f j+ 1 2 (µ ′ ) 1 ζ i+ 1 2 ln( f i+1 )(µ)-ln( f i )(µ) ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 ln( f j+1 )(µ ′ )-ln( f j )(µ ′ ) ∆ζ ′ j+ 1 2 (ln(ρ i+1 )-ln(ρ i )) ζ i+ 1 2 ∆ζ i+ 1 2 ζ 2 i+ 1 2 ζ ′2 j+ 1 2 ∆ζ i+ 1 2 ∆ζ ′ j+ 1 2 dµdµ ′ .
Next it remains to study the term

∑ i G 1 i+ 1 2 -G 1 i-1 2 ∆ζ i ∆ζ i α i = -∑ i G 1 i+ 1 2 (α i+1 -α i ) .
(5.12)

By using the discretized form of the collision operator Q 0 ( f 1 ) given by (5.4), the righthand side of (5.12) can be expressed as

-∑ i,j ζ i+ 1 2 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 )   F 0 j+ 1 2 1 ζ i+ 1 2 F 1 i+1 -F 1 i ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 F 1 i+ 1 2 F 0 j+1 -F 0 j ∆ζ ′ j+ 1 2   (α i+1 -α i )ζ ′2 j+ 1 2 ∆ζ ′ j+ 1 2 . (5.13)
The definitions of F 0

i+ 1 2 and F 1 i+ 1 2 lead to ∑ i G 1 i+ 1 2 -G 1 i-1 2 ∆ζ i ∆ζ i α i = -∑ i,j 1 -1 1 -1 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 )ζ 2 i+ 1 2 ζ ′2 j+ 1 2 f j+ 1 2 (µ ′ ) f i+ 1 2 (µ) ln( f i+1 )(µ)-ln( f i )(µ) ∆ζ i+ 1 2 ζ i+ 1 2 - ln( f j+1 )(µ ′ )-ln( f j )(µ ′ ) ∆ζ ′ j+ 1 2 ζ ′ j+ 1 2 1 ∆ζ i+ 1 2 ζ i+ 1 2 α i+1 µ-α i µ ∆ζ i+ 1 2 ∆ζ i+ 1 2 ∆ζ ′ j+ 1 2 dµdµ ′ .
Hence by using that ln( f i ) = ln(ρ i )-α i µ, we obtain

∑ i G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i ∆ζ i ln(ρ i )-∑ i G 1 i+ 1 2 -G 1 i-1 2 ∆ζ i ∆ζ i α i = -∑ i,j 1 -1 1 -1 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 ) f j+ 1 2 (µ ′ ) f i+ 1 2 (µ)ζ 2 i+ 1 2 ζ ′2 j+ 1 2   1 ζ i+ 1 2 ln( f i+1 )(µ)-ln( f i )(µ) ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 ln( f j+1 )(µ ′ )-ln( f j )(µ ′ ) ∆ζ i+ 1 2   1 ζ i+ 1 2 ln( f i+1 )(µ)-ln( f i )(µ) ∆ζ i+ 1 2 ∆ζ i+ 1 2 ∆ζ ′ j+ 1 2 dµdµ ′ .
By exchanging the indexes i and j and by using the symmetry of J, we get

∑ i G 0 i+ 1 2 -G 0 i-1 2 ∆ζ i ∆ζ i ln(ρ i )-∑ i G 1 i+ 1 2 -G 1 i-1 2 ∆ζ i ∆ζ i α i = - 1 2 ∑ i,j 1 -1 1 -1 J(ζ i+ 1 2 ,ζ ′ j+ 1 2 ) f j+ 1 2 (µ ′ ) f i+ 1 2 (µ)ζ 2 i+ 1 2 ζ ′2 j+ 1 2   1 ζ i+ 1 2 ln( f i+1 )(µ)-ln( f i )(µ) ∆ζ i+ 1 2 - 1 ζ ′ j+ 1 2 ln( f j+1 )(µ ′ )-ln( f j )(µ ′ ) ∆ζ ′ j+ 1 2   2 ∆ζ i+ 1 2 ∆ζ ′ j+ 1 2 dµdµ ′ .

Numerical tests

Now, we turn considering numerical applications to illustrate the properties of the previous method. For this numerical application, we chose a one dimensional direction of propagation. Therefore we take µ∈[-1,1] as the direction of propagation variable instead of Ω. Two main numerical test cases are presented in this section. The first one, developed in the two first subsections, is devoted to the the full discrete model

         1 ∆t f 0,N+1 i -f 0,N i = Q 0 i , 1 ∆t f 1,N+1 i -f 1,N i = Q 1 i , (6.1)
where the operators Q 0 i and Q 1 i are defined in (5.3). The scaling used here gives time steps ∆t of the same order of the electron-ion collision time.

In the first subsection we consider the convergence from an uncentered Maxwellian to a centered one. Next the second test case illustrates the convergence from a bi-Maxwellian distribution function to a centered one. The second aspect concerns the good working and the accuracy of the method completed with the transport part in one space dimension. We illustrate in the third subsection the propagation of an electron beam in vacuum which converges to a steady state.

Convergence to the steady state with an uncentered Maxwellian

This test case is devoted to study the evolution of an uncentered distribution function to the steady state thanks to the collision effects. The initial distribution function f is a Maxwellian distribution function centered around ζ d = 2 with a density n = 1 m -3 and a temperature T = 1 KeV. More precisely f is defined by

f (t = 0,ζ) = m 2π 3/2 exp - (ζ -ζ d ) 2 2 .
The energy domain is discretized with 64 points and the length of the energy space is denoted by ζ max =8 (for the 4 thermal velocity). We choose a time step satisfying the CFL condition. The scaling used here gives time steps of the order of the electron-ion collision time. The numerical results are computed at time t = 100. For the validation procedure, the numerical results are compared with the exact solution, i.e. the centered Maxwellian distribution with the same mass and energy than f (t = 0,ζ). The boundary conditions for the energy domain satisfy also naturally the positivity constraints of f 0 and the conservation of the mass, i.e. at the left boundary f 0 -1 = f 0 0 and f 1 -1 =f 1 0 .

Anisotropic case

In this subsection, an anisotropic initial function is chosen such as

       F 0 (t = 0,ζ) = exp(-(ζ -2) 2 ) , F 1 (t = 0,ζ) = F 0 (t = 0,ζ) 3 .
Figures 2 (a) and 2 (b) show the conservation of mass and energy in time which is consistent with the conservation laws given in subsection 3.1. Numerically the conservation is not perfect, because we considered a finite energy domain. Therefore the most energetic particles leave the domain during the simulation. Moreover, the evolution of the entropy induced by the scheme is decreasing in time (2 (c)) as shown by Theorem 4.1.

From figures 3 (a) and 3 (b) the realizability domain is still preserved during the simulation. Moreover the isotropic part F 0 converges to a centered Maxwellian. This fact can be explained by the electron-electron collision operator effect whereas the anisotropic part F 1 converges to null state because of the electron-ion collision operator effect. From figure 3 (c), the scheme has converged quickly from t = 80. Figure 3 (d) represents the evolution of the distribution function F with µ = 1, reconstructed from its isotropic part F 0 and its anisotropic part F 1 . We remark that the distribution function F converges to a centered Maxwellian.

In order to compare the M 1 model with the P 1 model ( [START_REF] Grad | On Kinetic Theory of the Rarefied Gases[END_REF]), consider the initial data [START_REF] Buet | Numerical analysis of conservative and entropy schemes for the Fokker-Planck-Landau equation[END_REF] shows that the M 1 model preserves realizability domain whereas the P 1 model does not conserve the realizability domain from first iterations. Now to evaluate the computational time of the M 1 model, we compare M 1 results with results obtained with the kinetic code KETS ( [START_REF] Duclous | Modélisation et Simulation Numérique multi-échelle du transport cinétique électronique[END_REF][START_REF] Duclous | High order resolution of of the Maxwell-Fokker-Planck-Landau model intended for ICF application[END_REF]) described in introduction. For the simulation with the kinetic model, each component of the velocity domain is discretized with 16 points. The numerical test is performed with v x,max = v y,max = v z,max = 8. M 1 model runs during 305 s whereas the CPU time describing KETS simulation is about 8717 s. Therefore, M 1 simulation runs 28 times faster than the KETS one.

F 1 (t = 0,ζ) = F 0 (t = 0,ζ)/2 instead of F 1 (t = 0,ζ) = F 0 (t = 0,ζ)/3. Figure
= ∑ i |F 0,N+1 i -F 0,N i | and d f 1 = ∑ i |F 1,N+1 i -F 1,N i | (
We have also checked for the isotropic case, i.e. F 1 (t = 0,ζ) = 0, that the conservation laws are preserved. Indeed, we observe that mass and energy are conserved. Besides, the evolution of the entropy induced by the scheme is decreasing in time. Moreover, by plotting the evolution moments F 0 and F 1 , we observe that F 1 (t,ζ) remains equal to 0 and the isotropic part F 0 (t,ζ) converges to the same centered Maxwellian than for the anisotropic case. Notice that the scheme has converged quickly from t = 40. We remark also that the distribution function F (defined with F 0 (t,ζ) and F 1 (t,ζ)) converges to a centered Maxwellian distribution which is close to the exact equilibrium state with an L 2 relative error about 10 -4 . This enhances the good accuracy of the scheme.

Convergence to the steady state with a bi-Maxwellian distribution function

The initial data is now a bi-Maxwellian distribution function i.e. a sum of two Maxwellian functions with a density n = 1 m -3 and T = 1 KeV. The first one is centered around ζ d1 = 2 and the second one around ζ d2 = 6. So the bi-Maxwellian function is defined as

           F 0 (t = 0,ζ) = m 2π 3/2 exp - (ζ -ζ d1 ) 2 2 +exp - (ζ -ζ d2 ) 2 2 , F 1 (t = 0,ζ) = F 0 (t = 0,ζ) 3 . (6.2)
For this test case, the energy domain is discretized with 64 points. The numerical test is performed with ζ max = 20 which represents the length of the energy space. The scaling used here gives time steps to be of the order of the electron-ion collision time as previously. We choose the time step in order to satisfy the CFL condition and the results are computed at time t = 100. The boundary conditions for the energy domain satisfies also naturally the positivity constraints of F 0 and the conservation of the mass, i.e. at the left boundary F 0 -1 = F 0 0 and F 1 -1 = -F 1 0 . Figures 5 (a Figures 6 (a) and 6 (b) confirm that the realizability domain is preserved. Moreover the isotropic part F 0 converges to a centered Maxwellian whereas the anisotropic part F 1 converges to null state. As shown on figure 6 (c), the scheme has converged quickly from In order to test the robustness of our scheme for the conservation of the realizability domain we consider initial conditions close to the boundary of the realizability domain as

           F 0 (t = 0,ζ) = m 2π 3/2 exp - (ζ -ζ d1 ) 2 2 +exp - (ζ -ζ d2 ) 2 2 , F 1 (t = 0,ζ) = 9F 0 (t = 0,ζ) 10 .
(6.3)

Figure 7 (a) shows that the realizability domain is still preserved by the numerical scheme. Moreover, from 7 (b), we observe that the approched solution relaxes toward the correct asymptotic steady state.

Evolution of an electron beam

Electron transport is an important process in plasma physics which impacts on many areas including laser-plasma. That is why now we turn considering the full discrete model including transport part and collision part to illustrate the good behavior and the accuracy of the method on a one space dimension domain. To discretize the transport equation, we use a Finite Volume method coming from the HLL method introduced by Harten and al. ( [24]). The full discrete model can be written such as

       1 ∆t f 0,N+1 i -f 0,N i + 1 ∆x F N i+ 1 2 -F N i-1 2 = Q 0 i , 1 ∆t f 1,N+1 i -f 1,N i + 1 ∆x G N i+ 1 2 -G N i-1 2 = Q 1 i ,
where F , G are HLL flux and the operators Q 0 j and Q 1 j are defined in (5.3). This test case consists to study the convergence of an electron beam, whose atomic mass is Z = 10, in the vacuum to a steady state. The initial distribution function in the one dimension space domain is null. The boundary condition on the right side is a Dirichlet condition chosen in order to satisfy realizability domain and to create an electron beam

F 0 -1 = m 2π 3 2 exp - (ζ -ζ d ) 2 2 and F 1 -1 = f 0 -1 3 .
For this test case, the energy domain is discretized with 64 points and the space domain is discretized with 100 points. The numerical test is performed with ζ max =8 and x max =5. The scaling used here gives time steps to be of the order of the electron-ion collision time.

We choose a time step satisfying the CFL condition. The numerical results are computed ). Then at time t = 0.15, figures 9 (c) and 9 (d) show that the isotropic part F 0 begins to relax to a Maxwell function centered around ζ = 0 while the anisotropic part F 1 is close to 0 at the end of the space domain. At the end of the computation, at time t=300, the isotropic part has converged to a Maxwellian function centered around 0. This Maxwellian is not uniform in space (figure 8). This phenomenon can be explained by the fact that the electron beam is still injected at the beginning of the space domain. Concerning the anisotropic part, F 1 is null except at the beginning of the space domain (8) for the same reason as for the isotropic part. The results obtained transcribe the convergence to the steady state due to the collision operators. Notice also that the realizability domain is still preserved in time.

= ∑ i |F 0,N+1 i -F 0,N i | and d f 1 = ∑ i |F 1,N+1 i -F 1,N i | (
In order to compare the entropic and the arithmetic averages, we consider the test case of the electron beam, with a less refined mesh. We take 16 points in the velocity variable for ζ max =8 and 30 points in the space variable. Figure 10 represents the absolute value of the anisotropy factor α= f 1 f 0 at time 40.4. We observe that the realizability domain is preserved when the entropic average is considered whereas it is not conserved for the arithmetic average. To evaluate the good behavior and the accuracy of the M 1 model, we compare M 1 results with results obtained with the kinetic code called KETS ( [START_REF] Duclous | Modélisation et Simulation Numérique multi-échelle du transport cinétique électronique[END_REF][START_REF] Duclous | High order resolution of of the Maxwell-Fokker-Planck-Landau model intended for ICF application[END_REF]) described in introduction. For the simulation with the kinetic model, each component of the velocity domain is discretized with 16 points and the space domain is discretized with 100 points. The numerical test is performed with v x,max = v y,max = v z,max = 8 and x max = 5. We choose a time step satisfying the CFL condition. The numerical results are computed at time t = 100. Figure 11 represents the mass with respect to the space domain computed with the M 1 model and with the KETS code. We observe that the two masses obtained with the two codes lead to constant distribution function at the end of the space computation domain. This fact is explained by the convergence of the distribution function toward a centered Maxwellian. Moreover it can be shown that the M 1 model creates a boundary layer of length 1 3 . That is why the two results are different on [0, 1 3 ]. On the remaining area, results are observed to be close. Moreover it has not been possible to compute the KETS code for a more refined mesh because of time computing.

Conclusion

We have developed a new model for electron transport in plasma taking into account the collisions between electrons and between electrons and ions. Its construction is essentially based on an entropy minimization principle and an integration of the kinetic model w.r.t the angle variable of the microscopic velocity. Next an entropic semi-discretization in energy of the model based on an entropic average is proposed. Finally numerical simulation in the context of the full discretized model are performed. The results are in a very good agreement with the properties proved on the semi-discrete model and results obtained with kinetic model. Several extension and application of this work can be considered. For instance, the generalization of the present results to a N moments model is postponed to a future paper.

A Computation of the B coefficient

The coefficient B (3.9) can be computed as

Figure 2 :

 2 Figure 2: Evolution of mass (a), energy (b) and entropy (in log scale) (c) in the case of ζ max = 8 with 64 points for anisotropic initial data.

Figure 3 :

 3 Figure 3: Description of the moments F 0 and F 1 with respect to the energy ζ at initial time (a) and at final time (b), evolution of the convergence in logarithmic scale i.e. d f 0 = ∑ i |F 0,N+1

  c) and evolution of the distribution function F with µ = 1 (d) in the case of ζ max = 8 with 64 points for anisotropic data.

Figure 4 :

 4 Figure 4: Description of the moments F 0 and F 1 with respect to the energy ζ at time t = 0.15 with ζ ∈ [0 : 8] for the M 1 model (a) and with ζ ∈ [0 : 8] for the P 1 model (b).

  ) and 5(b) show the mass and the energy conservation in time. The evolution of the entropy induced by the scheme is still decreasing in time as shown by Theorem 4.1.

Figure 5 :

 5 Figure 5: Evolution mass (a), energy (b) and entropy(c) in the case of ζ max = 20 with 64 points for bi-Maxwellian initial data.

Figure 6 :

 6 Figure 6: Description of the moments F 0 and F 1 with respect to the energy ζ at initial time (a) and at final time (b), evolution of the convergence in logarithmic scale i.e. d f 0 = ∑ i |F 0,N+1

  c) and evolution of the distribution function F in time with µ =1 (d) in the case of ζ max =8 with 64 points for bi-Maxwellian data.

Figure 7 :Figures 9

 79 Figure 7: Moments F 0 and F 1 with respect to the energy ζ at final time (a) and evolution of the distribution function F in time with µ = 1 (d) in the case of ζ max = 8 for the initial conditions 6.3

Figure 8 :Figure 9 :

 89 Figure 8: Description of moment F 0 (a) and F 1 (b) with respect to the energy ζ at time t =100 for x =0, x =0.5 and x = 5 in the case of ζ max = 8 with 64 points and x max = 5 with 100 points for electron beam data.

Figure 10 :

 10 Figure 10: Computation of |a|, a = f 1 / f 0 at time 40.4 with the arithmetic average (a) and the entropic average (b)

Figure 11 :

 11 Figure 11: Description of the mass ρ in space domain at final time obtained with the M 1 model and KETS model in the case of ζ max = 8 with 64 points and x max = 5 with 100 points for electron beam data.

  1 2 )∆ζ i . Let h i be an approximation of h(ζ i ) for all distribution function h and h i+ 1

	2	an approximation of h(ζ i+ 1 2

Acknowledgments

This work was partially supported by EURATOM within the "Keep-in-Touch"activities and was granted access to the HPC resources of CINES under the allocation 2011-056129 made by GENCI (Grand Equipement National de Calcul Intensif). We are thankfull with Vladimir Tikhonchuk for fruitfull discussions about the pertinence of the physical model.