Kim-Claire Le

Thanh -Septembre 2006 - Rapport Cea-R-6121 -Christophe Buet

Kim-Claire Le Thanh

Christophe Buet

ABOUT POSITIVE, ENERGY CONSERVATIVE AND EQUILIBRIUM STATE PRESERVING SCHEMES FOR THE ISOTROPIC

Sur les schémas positifs, conservant l'énergie et les états d'équilibre pour l'équation de Fokker-Planck-Landau isotrope» Résumé -Dans ce rapport on s'intéresse à la discrétisation des termes de collisions de l'équation de Fokker-Planck-Landau. En particulier on étudie les collisions entre des électrons isotropisés par un bain de particules lourdes, par exemple des ions. La discussion porte sur les schémas positifs qui conservent la masse, l'énergie et les états d'équilibre maxwelliens. En premier lieu on analyse en détails le schéma de Chang et Cooper pour ce terme de collisions non linéaires : construction du schéma, positivité et propriétés de conservation. On montre que certaines variantes de ce schéma, dérivées de l'écriture de l'opérateur de Fokker-Planck sous la forme Rosenbluth, peuvent ne pas être positives ou ne pas conserver l'énergie. Nous présentons aussi une nouvelle version du schéma de Chang et Cooper, écrite à partir de la forme Landau, qui possède toutes les bonnes propriétés. Pour terminer, on propose deux autres nouveaux schémas, eux aussi dérivés de la forme Landau, qui montrent que le schéma de Chang et Cooper n'est pas le seul à posséder de bonnes propriétés. Pour tous ces schémas, nous analysons clairement les propriétés de conservation de la densité et de l'énergie, la positivité des solutions et la conservation des états d'équilibre. On s'intéresse également aux cas où les collisions sont maxwelliennes ou coulombiennes.

ABOUT POSITIVE, ENERGY CONSERVATIVE AND EQUILIBRIUM STATE PRESERVING SCHEMES FOR THE ISOTROPIC FOKKER-PLANCK-LANDAU EQUATION par Christophe BUET et Kim-Claire LE THANH C E A / D A M I L E -D E -F R A N C E D É P A R T E M E N T S C I E N C E S D E L A S I M U L A T I O N E T D E L ' I N F O R M A T I O N S E R V I C E N U M É R I Q U E E N V I R O N N E M E N T E T C O N S T A N T E S D É P A R T E M E N T D E P H Y S I Q U E T H É O R I Q U E E T A P P L I Q U É E S E R V I C E P H Y S I Q U E D E S P L A S M A E T E L E C T R O M

The Fokker-Planck-Landau equation is commonly used in plasma physics when studying kinetic effects between charged particles under coulomb interaction. The isotropic Fokker-Planck-Landau operator is generally used in the modelling of inertial controlled fusion. More precisely, it is used to describe the electronic energy transport phenomena in laser produced plasma. In some conditions, it is well known that the fluid theory, for which the hydrodynamic equations are closed using the law for the thermal fluxes proposed by Spitzer-Harm [START_REF] Spitzer | Transport phenomena in a completely ionized gaz[END_REF], is not valid [START_REF] Epperlein | Fokker-Planck modelling of electrons transport in laser-produced plasmas[END_REF][START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF]. A more accurate solution is to use a model based on the expansion of the FPL operator in spherical harmonics and to retain the two first terms, and the isotropic FPL operator is the leading order term [START_REF] Epperlein | Fokker-Planck modelling of electrons transport in laser-produced plasmas[END_REF][START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF][START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF]. There are also applications in the astrophysical field where the FPLE is used for star's clusters modelling [START_REF] Cohn | Numerical integration of the Fokker-Planck equation and the evolution of stars clusters[END_REF][START_REF] Cohn | Late core collapse in star clusters and the gravothermal instability[END_REF]. The most popular scheme for Fokker-Planck type equations is the Chang and Cooper method [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF]. This method was originally devoted to linear Fokker-Planck equations and it was shown in [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF] that in this case this method is positive and preserves the equilibrium states. On the contrary a paper of Larsen et al. [START_REF] Larsen | Discretization methods for one-dimensional Fokker-Planck operators[END_REF] shows that this method applying for nonlinear Fokker-Planck equations could produce non positive solutions. For all that, this method is used for the isotropic Landau equation by Langdon [START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF], in the SPARK code by Epperlein [START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF], by Kingham and Bell [START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF] in the code IMPACT or by Cohn [11] in the astrophysical context. To our knowledge, there is no rigorous proof of the positivity and the conservation of the energy of this scheme when it is applied to isotropic Landau equation. For the isotropic Landau equation the Chang and Cooper scheme is derived from the Rosenbluth form of the equation which is not the most powerful form to check energy conservation and could lead to mistake at boundary. There exist also conservative, positive and entropy schemes for this equation [START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF][START_REF] Buet | Conservative and entropy schemes for the isotropic Fokker-Planck-Landau equation[END_REF][START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF][START_REF] Degond | An entropy scheme for the Fokker-Planck collision operator of plasma kinetic theory[END_REF][START_REF] Lemou | Multipole expansions for the Fokker-Planck-Landau operator[END_REF][START_REF] Lemou | Numerical Algorithms for Axisymmetric Fokker-Planck-Landau Operators[END_REF] based on the "Log" weak symmetric form of the Landau equation, but these methods suffer from some limitations due to the use of Logarithm of the distribution function or, when not using "Log", they can only be defined for uniform mesh in the velocity space.

Thus it would be interesting to clarify the situation about positive, conservative and equilibrium states preserving velocity discretization of the isotropic Landau equation. This is the aim of this work. In the most important part of this paper we focus our attention on the Chang and Cooper method. By starting from the weak symmetrized form of the isotropic Landau operator we construct two variants of the Chang and Cooper method. A new one, called S 1 , and an other one, called S 2 , that is described by Langdon et al. [START_REF] Decoster | [END_REF][START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF] in the Coulombian case but constructed directly from the Rosenbluth form of the operator. These two schemes only differ in the boundary condition. We prove that these two variants conserve the energy but we show that the variant which fits with the Langdon's scheme is not positive. Moreover we make also some remarks about simplifications of these schemes which could lead to the loss of the positivity or the loss of the energy conservation. These remarks would be useful for the time discretization, especially for time implicit discretization. We clarify also the boundary conditions at the end of the domain. As we will see this leads to the conclusion that, due to bad boundary conditions, the scheme used by Epperlein [START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF] or by Kingham and Bell [START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF], and also obtained using the Rosenbluth form, is not conservative in energy (because in collision operator we can't reverse the role of test particles and target particles) but it is positive.

In other hand we also show that the Chang and Cooper method is not the only way to obtain positive, conservative and equilibrium preserving states even on non-uniform meshes. We propose two new other schemes that share the same properties with the Chang and Cooper method. One of them, called S 3 , is based on the work of Larsen et al. [START_REF] Larsen | Discretization methods for one-dimensional Fokker-Planck operators[END_REF]. The other, called S 4 , is based on the "Log" form of the equation (it is a new version of the scheme described in [START_REF] Degond | An entropy scheme for the Fokker-Planck collision operator of plasma kinetic theory[END_REF] applied to the isotropic FPL equation). Thus it is also an entropy decaying scheme. In the case of uniform meshes these two schemes are nothing but the conservative positive and entropy scheme of Berezin et al. [START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF] and studied in detail by one of the authors [START_REF] Buet | Conservative and entropy schemes for the isotropic Fokker-Planck-Landau equation[END_REF][START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF]. This paper is organized as follow: in the first part, we recall the continuous FPL equation in the homogeneous and isotropic case and its properties. We recall also the different forms of the equation in the case of Maxwellian or Coulombian potentials. In a second part we present the two Chang and Cooper type schemes and the two non Chang and Cooper ones. We clearly show all the properties of these schemes: energy conservation, equilibrium states and positivity. We also indicate for each of them, the simplifications in the case of a uniform mesh, or in the maxwellian case. Some comments are also made about the advantages or disadvantages of each of these four schemes especially for the implementation.

Fokker-Planck-Landau equation

We present the homogeneous non-linear Fokker-Planck-Landau equation (FPL equation) in the isotropic case where the distribution function f ( x, v, t) depends only on the modulus of the velocity v = v and on the time t, in other words f ( x, v, t) = f (v, t). We consider f as a function of ε = v 2 where ε is the energy variable. In this case, Fokker-Planck-Landau equation can be written:

∂f (ε, t) ∂t = Q(f )(ε) = 1 √ ε ∂ ∂ε ∞ 0 g(ε, ε ) f (ε ) ∂f (ε) ∂ε -f (ε) ∂f (ε ) ∂ε dε , (2.1)
where g(ε, ε ) is positive, symmetric and increasing (g(ε, ε ) = min(ε

3 2 , ε 3
2 ) for Coulombian interactions and g(ε, ε ) = ε 3 2 ε 3 2 for Maxwellian interactions) and Q(f ) is the so-called FPL collision operator. This operator can be written in the following weak form (let φ(ε) be any test (smooth and decaying) time independent function)

∞ 0 ∂f (ε, t) ∂t φ(ε) √ εdε = ∞ 0 φ(ε) ∂ ∂ε ∞ 0 g(ε, ε )(f (ε ) ∂f (ε) ∂ε -f (ε) ∂f (ε ) ∂ε ) dε dε , (2.2
) and by integrating (2.2) by parts

∞ 0 ∂f (ε, t) ∂t φ(ε) √ εdε = - 1 2 ∞ 0 ∞ 0 ( ∂φ(ε) ∂ε - ∂φ(ε ) ∂ε ) g(ε, ε ) (f (ε ) ∂f (ε) ∂ε -f (ε) ∂f (ε ) ∂ε ) dε dε .
(2.3) Note that the FPL equation can be equivalently written in the so-called "Log" weak form

∞ 0 ∂f (ε, t) ∂t φ(ε) √ εdε = - 1 2 ∞ 0 ∞ 0 ( ∂φ(ε) ∂ε - ∂φ(ε ) ∂ε )g(ε, ε )f (ε )f (ε)( ∂ log f (ε) ∂ε - ∂ log f (ε ) ∂ε ) dε dε .
(2.4)

Properties

Let us recall the most important properties of the problem (2.1)

1. This operator satisfies the conservation of mass (respectively energy) by choosing φ

(ε) = 1 (respectively φ(ε) = ε) in (2.3) ρ = ∞ 0 f (ε, t) √ εdε , (2.5
)

ρE = ∞ 0 f (ε, t)ε 3 2 dε . (2.6)
Note that the conservation properties are a consequence of the symmetry property (between ε and ε ) of the collision operator. Let us also mention that the temperature T of the plasma is defined as 3 2 ρT = ρE.

2. Any function of the type ψ(ε) = α exp(-βε) where α and β are arbitrary constants (β > 0) is a stationary solution of equation (2.1). The laws of conservation select, for this two-parameter set of functions, the unique equilibrium solution corresponding to the initial condition f (ε, 0) = f 0 (ε).

The entropy defined by

H(f ) = ∞ 0 f (ε) log(f (ε)) √ ε dε , (2.7)
satisfies the classical H-Theorem. That means the entropy decays with time (by letting φ(ε) = log(f (ε)) in the weak formulation (2.4)) and

dH(f ) dt = 0 ⇐⇒ f (ε) = α exp(-βε) . (2.8)
This is formally equivalent to say that f is an equilibrium function, that is

Q(f ) = 0.
The FPL operator can be rewritten in the following diffusive form

Q(f )(ε) = 1 √ ε ∂ ∂ε (E(f ) f (ε) + D(f ) ∂f ∂ε ) , (2.9)
where

D(f ) = ∞ 0 g(ε, ε )f (ε )dε , and
E(f ) = - ∞ 0 g(ε, ε ) ∂f (ε ) ∂ε dε .
(2.10)

Remark 1. To preserve the number of particles and the energy, Bobylev and Chuyanov [START_REF] Bobylev | On the numerical solution of Landau's kinetic equation[END_REF] write the FPL operator in the form

Q(f )(ε) = 1 √ ε ∂ 2 ∂ε 2 W (f, ε) ,
where

W (f, ε) = ε 0 ∞ 0 g(ε , ε ) f (ε ) ∂f (ε ) ∂ε -f (ε ) ∂f (ε ) ∂ε dε dε , (2.11)
or, from (2.10)

W (f, ε) = ε 0 E(f ) f (ε ) + D(f ) ∂f (ε ) ∂ε dε .
(2.12)

By symmetry (between ε and ε in (2.11)) W (f, ε) vanishes at infinity and also at zero point, therefore density and energy are conserved.

Remark 2. The mathematical analysis of the Landau equation is done by Desvillettes and

Villani in [START_REF] Desvillettes | On the spatially homogeneous Landau equation for hard potentials. Part I: Existence, uniqueness and smoothness[END_REF][START_REF] Desvillettes | On the spatially homogeneous Landau equation for hard potentials. Part II: H-Theorem and applications[END_REF]. Concerning the derivation from the Boltzmann equation in the case of Coulombian potential there is the work of Degond and Lucquin [START_REF] Degond | The Fokker-Planck asymptotics of the Boltzmann collision operator in the Coulomb case[END_REF]. There is also the work of Desvillettes for the derivation from the Boltzmann equation when the collisions become grazing [START_REF] Desvillettes | On Asymptotics of the Boltzmann Equation when the Collisions Become Grazing[END_REF]. For the development of the FPL operator in spherical harmonics we refer to the work of Shkarofsky et al. [START_REF] Shkarofsky | The Particle Kinetics of Plasmas[END_REF].

Coulombian interactions

In the Coulombian case the collisions coefficients E(f ) and D(f ) describe drift and diffusion between electrons and are given by

E(f ) = - ε 0 ε 3 2 ∂f (ε ) ∂ε dε -ε 3 2 ∞ ε ∂f (ε ) ∂ε dε = - ε 0 ε 3 2 ∂f (ε ) ∂ε dε + ε and D(f ) = ε 0 ε 3 2 f (ε ) dε + ε 3 2 ∞ ε f (ε ) dε , (2.14)
respectively. Integrating by parts and using the fact that f vanishes at ∞ we obtain

E(f ) = 3 2 ε 0 √ ε f (ε ) dε , (2.15)
and

D(f ) = 3 2 ε 0 √ ε ∞ ε f (ε ) dε dε . (2.16)
Deriving once E(f ) and twice D(f ) we have

f (ε) = 2 3 1 √ ε ∂E(f ) ∂ε , (2.17)
and

f (ε) = - 2 3 ∂ ∂ε 1 √ ε ∂ ∂ε D(f ) . ( 2
W (f, ε) = ε 0 3 2 ε 0 √ ε f (ε )dε f (ε ) dε + ε 0 D(f, ε ) ∂f (ε ) ∂ε dε .
By integrating by parts, we can write

W (f, ε) = 3 2 ε 0 √ ε f (ε ) ∞ ε f (ε )dε dε - 3 2 ε 0 √ ε f (ε )dε ∞ ε f (ε ) dε - ∞ 0 ∂D(f, ε ) ∂ε f (ε )dε + D(f, ε)f (ε) . As ∂D(f, ε) ∂ε = 3 2 √ ε ∞ ε f (ε )
dε the first and the third integral in the right-hand side of the previous equation become identical, thus we obtain

W (f, ε) = - 3 2 ε 0 √ ε f (ε ) dε ∞ ε f (ε ) dε + D(f, ε) f (ε) , in other words W (f, ε) = 2 3 - E(f, ε) √ ε ∂ ∂ε D(f, ε) + D(f, ε) √ ε ∂ ∂ε E(f, ε) .

Maxwellian interactions

In the Maxwellian interactions case we simply get

E(f ) = 3 2 ε 3 2 ρ and D(f ) = ε 3 2 ρE = 3 2 ε 3 2 ρ T .
Then, equation (2.9) is linear (in the sense that E(f ) and D(f ) are independent of f ). We can write

Q(f )(ε) = 3 2 ρ √ ε ∂ ∂ε ε 3 2 ( f (ε) + T ∂f (ε) ∂ε ) .
3 Semi-discretized problem

In this section we focus on the discretization in the energy variable. For numerical simulations, we reduce the integration domain in FPL equation to a bounded domain in the variable ε where ε ∈ [0, E]. Thus we consider the approximate problem of (2.3) defined by

E 0 ∂f (ε, t) ∂t φ(ε) √ εdε = - 1 2 E 0 E 0 ( ∂φ(ε) ∂ε - ∂φ(ε ) ∂ε ) g(ε, ε ) (f (ε ) ∂f (ε) ∂ε -f (ε) ∂f (ε ) ∂ε ) dε dε . (3.1) Let us introduce {ε i } 1≤i≤N +1 an increasing sequence such that ε 1 = 0, ε N = E and ∆ε i+ 1 2 = ε i+1 -ε i .
We suppose that {∆ε i+ 1 2 } 1≤i≤N is constant or increasing sequence. We define v i+ 1 2 as the mean value of the velocity on [ε i , ε i+1 ] i.e.

v i+ 1 2 = 1 ∆ε i+ 1 2 ε i+1 ε i √ ε dε = 2 3∆ε i+ 1 2 (ε 3 2 i+1 -ε 3 2 i ) , so, we get ε i+ 1 2 and can define (v 2 ∆v) i = 1 2 ε i+ 1 2 ε i-1 2 √ εdε = 1 2 √ ε i ∆ε i = 1 3 (ε 3 2 i+ 1 2 -ε 3 2 i-1 2 ) . (3.2)
By convention we set ε 1 2 = 0. Any function f (ε, t) is approximated on the grid by values {f i } 1≤i≤N supposed be approximations of {f (ε i )} 1≤i≤N . We also introduce the notations

(∆φ) i+ 1 2 = φ i+1 -φ i and (Dφ) i+ 1 2 = (∆φ) i+ 1 2 (∆ε) i+ 1 2
as an approximation of the partial derivative

∂ ε φ(ε i+ 1 2
).

First, we consider the discretization of the expression

E 0 ∂f ∂t φ(ε) √ εdε for any function φ.
By writing

E 0 ∂f ∂t φ √ εdε N i=1 ε i+ 1 2 ε i-1 2 ∂f ∂t φ(ε) √ ε dε,
and using the standard quadrature formula with respect to the measure √ εdε, we approximate it by

N i=1 (φ i df i dt ε i+ 1 2 ε i-1 2 √ ε dε) def = N i=1 c i φ i df i dt (3.3)
Thus the weights c i are defined by

c 1 = 2 3 ε 3 2 3 2 , c i = 2 3 (ε 3 2 i+ 1 2 -ε 3 2 i- 1
2

) for i = 2, • • • , N .
Below, we present various strategies to construct schemes that have properties of conservation, and if possible positivity and entropy decaying. These schemes differ intrinsically in the way we discretize the right-hand side of (3.1):

(r.h.s.) = -1 2

N -1 i=1 N -1 j=1 ε i+1 ε i ε j+1 ε j ( ∂φ(ε) ∂ε - ∂φ(ε ) ∂ε ) g(ε, ε ) (f (ε ) ∂f (ε) ∂ε -f (ε) ∂f (ε ) ∂ε ) dε dε .
(3.4) We must also introduce the Maxwellian associated to a distribution function f . For a distribution function f we define ρ and ρE the discretized analogous of density (2.5) and energy (2.6) as

ρ = N i=1 c i f i , ρE = N i=1 ε i c i f i .
The temperature is still defined by 3 2 ρT = ρE. For a distribution function f , we denote by M = α exp(-βε) the Maxwellian which has the same mass and energy as f . It's easy to check that

T = N i=1 ε i c i M i N i=1 c i M i = N i=1 ε i c i exp(-βε i ) N i=1 c i exp(-βε i ) is a strictly monotone (decreasing) function of β with lim β→0 T = T max = N i=1 ε i c i N i=1 c i and lim β→+∞ T = T min = 0.
Thus for any distribution function f such that 0 ≤ T ≤ T max there is a unique β ≥ 0 such that T = T and consequently M is unique. Note that β 1 T . For the rest of this work we consider only distribution functions such that 0 ≤ T ≤ T max , that is, we exclude distribution function for which M is an increasing function of ε (β ≤ 0)( that means we exclude distribu-

tion function such that T max ≤ T ≤ 2 3 E).
To check the positivity of the schemes we will have the need of the following well-known result:

Lemma 1. Consider the Cauchy problem for the ordinary differential equation

df dt = Lf f (t = 0) = f 0 with f = {f i } 1≤i≤N and with the square matrix L = L(f, t) such that L ij ≥ 0 for i = j, L ii ≤ 0.
If there exists a constant C such that ∀ i, j |L ij (f, t)| ≤ C then, if the solution exists on a time interval [0, T ] and if the initial data is positive i.e. f 0 > 0, the solution f is positive i.e. f i (t) > 0 for all i and for all time t ∈ [0, T ].

If there exists a constant C such that ∀ i, j |L ij (f, t)| ≤ C then, if the solution exists on a time interval [0, T ] and if the initial data is non-negative i.e. f 0 ≥ 0, the solution f is non-negative i.e. f i (t) ≥ 0 for all i and for all time t ∈ [0, T ].

Such an assumption for the matrix L ensures the existence and the uniqueness of a global solution in time if the matrix L is locally Lipchitz in f . If the matrix L is only continuous in f we have also the existence of a global solution in time but not the uniqueness. One can see [START_REF] Crouseix | Analyse numérique des équations différentielles[END_REF] for these well known results. Let us also recall that the solution in these cases is C 1 in the time variable.

Proof. First we suppose that f 0 > 0. Let t 0 ≤ τ be the first time for which there exists an index i 0 such that f i 0 (t 0 ) = 0. In [0, t 0 ] for all i, we have df i dt ≥ -Cf i and by integrating

f i (t) ≥ f 0 i exp(-Ct).
Thus the solution cannot vanish in finite time. The extension of this result to the case of a non-negative initial data (f 0 i ≥ 0 for all i) depends on the regularity of the matrix L with respect to f . If L is Lipschitz in f the solution depends continuously on the initial data, see by example [START_REF] Crouseix | Analyse numérique des équations différentielles[END_REF], and thus by taking a sequence f 0 ν > 0 such that f 0 ν → f 0 and by using the result above concerning strictly positive initial data, f is non-negative.

We suppose now that L is just continuous in f . First we assume that for all f and for all i, |L ii | > j =i L ij that is L(f ) is strictly diagonally dominant. We suppose that the solution becomes negative for a set I of indexes up to a time t 1 > 0 (I ⊂ {i such that f 0 i = 0}). Since the solution is continuous in time, min i f i is defined and continuous in time in [0, t 1 ]. Let J ⊂ I such that f j = min i f i for all j ∈ J and for all t ∈ [0, t 2 ], t 2 ≤ t 1 . We have the following lowerbound for j ∈ J:

f j (t 1 ) = t 1 0 ( k =j L jk (f )f k + L jj (f )f j )dt ≥ t 1 0 ( k =j,k∈I L jk (f )f k + L jj (f )f j )dt > 0
thus this is in contradiction with our assumption that in [0, t 1 ], f j < 0. Thus the solution verifies f ≥ 0. Now if L(f ) is not strictly diagonally dominant but verifies |L ij (f, t)| ≤ C, we can choose a constant K such that L(f ) -KId is strictly diagonally dominant. We have

df dt = Lf ⇐⇒ df dt -Kf = (L -KId)f,
where Id is the identity matrix. Thus if we set g = exp(-Kt)f , g verifies df dt = (L(g exp(Kt)) -KId)g = L(g)g,

and the matrix L(g) is strictly diagonally dominant. Using the above result for strictly diagonally dominant matrix we have g ≥ 0 and thus f ≥ 0. If we do not suppose that L(f ) is continuous in f , but if we suppose the existence of a solution continuous in time in a time interval [0, T ], the proof above is still valid on this time interval. Thus the solution is non-negative on [0, T ]. This ends the proof.

Chang and Cooper type schemes

One of the most popular method used for Fokker-Planck equations is due to Chang and Cooper [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF]. Originally this method was proposed for the linear Fokker-Planck equation and the construction is entirely devoted to the preservation of equilibrium states. The authors also show that in the linear case (and only in this case) the method provides also non-negative solutions. At contrary others authors, see [START_REF] Larsen | Discretization methods for one-dimensional Fokker-Planck operators[END_REF], have shown that for the non-linear Fokker-Planck equation this method can produce non-positive solutions. Let us recall the spirit of the Chang and Cooper method on the simple Fokker-Planck equation for a distribution function

f (v, t), v ∈ R, ∂f ∂t = ∂F ∂v = ∂ ∂v (vf + σ ∂f ∂v ),
that can be put under the form

∂f ∂t = σ ∂ ∂v (M ∂ ∂v ( f M )),
where M is a Maxwellian, and M = M (v) = exp(-|v| 2 /2σ). Note that M is the equilibrium, the long time behavior of the distribution function. When f = M , F vanishes. On a uniform grid of velocity space step ∆v and v i = i∆v, the Chang and Cooper method consists to discretize the diffusion as usual and the drift in such a way that for f = M the fluxes F are all equal to zero. One takes at cell interface v i+ 1 2

= (i + 1 2 )∆v F i+ 1 2 = σ ∆v (f i+1 -f i ) + v i+ 1 2 (1 -δ i+ 1 2 )f i+1 + δ i+ 1 2 f i ,
and forcing the fluxes to be zero for equilibrium leads to the following definition of δ i+ 1 2

δ i+ 1 2 = σ v i+ 1 2 ∆v - 1 exp( v i+ 1 2 ∆v σ ) -1 .
The scheme writes as

df i dt = 1 ∆v (F i+ 1 2 -F i-1 2
).

As explained in [START_REF] Buet | A Conservative and Entropy Scheme for a Simplified Model of Granular Media[END_REF][START_REF] Buet | About the Chang and Cooper method for linear Fokker-Planck equations[END_REF] this scheme is in fact an entropic scheme for the above example.

This method is also one of the most used for FPL equation. And this is a non-linear problem. None of the main work in this area [START_REF] Cohn | Numerical integration of the Fokker-Planck equation and the evolution of stars clusters[END_REF][START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF][START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF][START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF] contains the proof of the positivity of the Chang and Cooper method or the proof of the energy conservation. Moreover it is not clear that equilibrium states are preserved by this method for the FPL equation since the definition of the coefficients of Chang and Cooper, see [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF] of these coefficients, are defined in an implicit manner.

In this section we propose a new derivation of this method which intrinsically contains the conservation of the energy. We recall also the scheme developed by Langdon et al. [START_REF] Decoster | [END_REF][START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF], but in the case where potentials are Coulombian, and extend it to the general case. We show clearly why our new method provides non-negative solution and why equilibrium states are still preserved. As we see later, this analysis clarifies the confusion about the boundary condition used at one end of the domain of computation, more precisely in ε = E. In particular this analysis shows that, due to boundary condition at ε = E, we can have some doubt about the energy conservation for the version of the scheme used by Epperlein in [START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF] and also used by other people, [START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF] for example. The discretization is only in energy, but the analysis provides also some results about implicit time discretization, more precisely about some simplifications of the totally implicit time discretization which could be made: the only positive, conservative and equilibrium state preserving time implicit discretization of the Chang and Cooper method for the FPL equation is the fully implicit one.

As example, the schemes are applied to the case where potentials are Coulombian or Maxwellian.

A new variant: scheme S 1

Using for each integrals of (3.4) a midpoint quadrature formula, we approximate it by (r.h.s.) = -1 2

N -1 i=1 N -1 j=1 ((Dφ) i+ 1 2 -(Dφ) j+ 1 2 ) g i+ 1 2 ,j+ 1 2 (f j+ 1 2 (Df ) i+ 1 2 -f i+ 1 2 (Df ) j+ 1 2 ) ∆ε i+ 1 2 ∆ε j+ 1 2 , (3.5) or (r.h.s.) = - N -1 i=1 (Dφ) i+ 1 2 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 (f j+ 1 2 (Df ) i+ 1 2 -f i+ 1 2 (Df ) j+ 1 2 ) ∆ε j+ 1 2 , with g i+ 1 2 ,j+ 1 2 = g(ε i+ 1 2 , ε j+ 1 2
). Hence, the weak formulation of the semi-discretized model reads

N i=1 c i ∂f i ∂t φ i = - N -1 i=1 (Dφ) i+ 1 2 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 (f j+ 1 2 (Df ) i+ 1 2 -f i+ 1 2 (Df ) j+ 1 2 ) ∆ε j+ 1 2 . (3.6)
To simplify we note

K i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 (f j+ 1 2 (Df ) i+ 1 2 -f i+ 1 2 (Df ) j+ 1 2 ) ∆ε j+ 1 2
the numerical flux. By factorizing the term φ i in the right-hand side of (3.6), we have

N i=1 c i ∂f i ∂t φ i = N -1 j=2 (K j+ 1 2 -K j-1 2 ) φ j + φ 1 K 3 2 -φ N K N -1 2 .
Finally, by identifying the terms involving φ i in (3.6), we obtain the system of ordinary differential equation

df i dt = Q S 1 i 1 ≤ i ≤ N (3.7) with Q S 1 1 = K 3 2 /c 1 , Q S 1 i = (K i+ 1 2 -K i-1 2 )/c i for 2 ≤ i ≤ N -1 and Q S 1 N = -K N -1 2
/c N . We can rewrite the numerical flux

K i+ 1 2 = -f i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 + (Df ) i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 , 1 ≤ i ≤ N -1 (3.8)
and recognize the discretized collisions terms (2.10)

E i+ 1 2 = - N -1 j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 and D i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 , (3.9) therefore K i+ 1 2 = E i+ 1 2 f i+ 1 2 + D i+ 1 2 f i+1 -f i ∆ε i+ 1 2 , 1 ≤ i ≤ N -1 (3.10)
Remark 6. By integrating by parts the drift term reads

E i+ 1 2 = N -1 j=1 (g i+ 1 2 ,j+ 1 2 -g i+ 1 2 ,j-1 2 ) f j -g i+ 1 2 ,N -1 2 f N .
We assume that {g i+ 1 2 ,j+ 1 2 } 1≤i≤N is an increasing sequence. Then, even if the f i 's are positive,

E i+ 1 2 can be negative if f N = 0. In other hand, since the f i+ 1 2 's are positive, D i+ 1 2 remains positive.
We turn now to the discretization of

f i+ 1 2 . We suppose that f i+ 1 2 is an approximation of f (ε i+ 1 2
) given by the following definition [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF] Definition 1. The Chang and Cooper average f i+ 1 2 of quantities f i and f i+1 is defined by

f i+ 1 2 = δ i+ 1 2 f i + (1 -δ i+ 1 2 ) f i+1 1 ≤ i ≤ N -1 , (3.11) with δ i+ 1 2 = 1 α i+ 1 2 - 1 exp(α i+ 1 2 ) -1 1 ≤ i ≤ N -1 , (3.12)
and

α i+ 1 2 = E i+ 1 2 D i+ 1 2 ∆ε i+ 1 2 .
To be comfortable we denote h

(α) = 1 α - 1 (exp(α) -1) (thus δ i+ 1 2 = h(α i+ 1 2
)). This smooth function is decreasing and bounded such that h(-∞) = 1, h(0) = 1/2 and h(+∞) = 0. Moreover, his derivative is negative and bounded too (h (±∞) = 0). Note that whatever the value of α i+ . We recall that f is the N-dimensional column vector with components

{f i } 1≤i≤N . First α i+ 1 2 (for 1 ≤ i ≤ N -1) reads α i+ 1 2 = E i+ 1 2 D i+ 1 2 ∆ε i+ 1 2 =        - N -1 j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2        ∆ε i+ 1 2 , α i+ 1 2 =        - N -1 j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 ((1 -δ j+ 1 2 )f j+1 + δ j+ 1 2 f j ) ∆ε j+ 1 2        ∆ε i+ 1 2 ,
and

δ i+ 1 2 = h(α i+ 1 2
). Thus, assuming that δ is the (N-1)-dimensional column vector with com-

ponents {δ i+ 1 2 } 1≤i≤N -1 , we can write δ i+ 1 2 as a function of f and introducing {H i+ 1 2 } 1≤i≤N -1 we get δ i+ 1 2 (f ) = H i+ 1 2 (δ 3 2 , ..., δ k+ 1 2 , ..., δ N -1 2 , f 1 , ..., f k , ... f N ) ,
in other words, the system to be solved is

δ(f ) = H(δ(f ), f ) , (3.13)
where f is a distribution function that has density ρ and energy ρE and H is a (N-1)dimensional vector valued function. Before solving this non-linear problem, we have to prove the existence of solution. To begin with, as f is a distribution function that has density ρ and energy ρE,

{E i+ 1 2 } 1≤i≤N -1 and {D i+ 1 2
} 1≤i≤N -1 are well defined and continuous in respect to δ. Hence {α i+ 1 2 } 1≤i≤N -1 exists (positive or negative) and consequently H lies in

[0, 1] N -1 . Because h is continuous, H is a continuous mapping which carries [0, 1] N -1 into
itself and thanks to Brouwer's theorem H has a fixed point (i.e. there exists an δ(f ) with

H(δ(f ), f ) = δ(f )).
Unfortunately we don't get actually an effective proof for the uniqueness of solution (as, for example, the Picard's theorem). Nevertheless, we assume that we have sufficient conditions for the uniqueness of solution of the non-linear operator and that under these conditions the sequence generated by Newton's method converges to the solution. In each step of Newton's process we have to solve a linear system which associated matrix is non-sparse. That requires O(N 2 ) operations. But we will see that for Maxwellian and Coulombian potential this cost can be reduce to O(N ) operations.

Proposition 1. The flux K i+ 1 2 satisfies the following relation

K i+ 1 2 = A i+ 1 2 f i+1 -B i+ 1 2 f i , ∀i ; 1 ≤ i ≤ N -1 (3.14)
where

A i+ 1 2 = v(α i+ 1 2 ) D i+ 1 2 ∆ε i+ 1 2 , (3.15)
and

B i+ 1 2 = u(α i+ 1 2 ) D i+ 1 2 ∆ε i+ 1 2 , (3.16
)

with u(α) = α exp(α) -1 and v(α) = α exp(α) exp(α) -1 .
Proof. The proof is the same as in Chang and Cooper's paper [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF] but in the non-linear case. By substituting the Chang and Cooper average f i+ 1 2 in the right-hand side of the relation (3.10) we get

K i+ 1 2 = E i+ 1 2 ((1 -δ i+ 1 2 ) f i+1 + δ i+ 1 2 f i ) + D i+ 1 2 f i+1 -f i ∆ε i+ 1 2 ∀i ; 1 ≤ i ≤ N -1 . As α i+ 1 2 = E i+ 1 2 D i+ 1 2 ∆ε i+ 1 2 , we can write E i+ 1 2 = α i+ 1 2 D i+ 1 2 ∆ε i+ 1 2
. Thus

K i+ 1 2 = α i+ 1 2 D i+ 1 2 ∆ε i+ 1 2 ((1 -δ i+ 1 2 )f i+1 + δ i+ 1 2 f i ) + D i+ 1 2 f i+1 -f i ∆ε i+ 1 2
, that means

K i+ 1 2 = (α i+ 1 2 (1 -δ i+ 1 2 ) + 1) D i+ 1 2 ∆ε i+ 1 2 f i+1 + (α i+ 1 2 δ i+ 1 2 -1) D i+ 1 2 ∆ε i+ 1 2 f i .
Finally, if we develop δ i+ 1 2

we have

K i+ 1 2 = α i+ 1 2 exp(α i+ 1 2 ) exp(α i+ 1 2 ) -1 D i+ 1 2 ∆ε i+ 1 2 f i+1 - α i+ 1 2 exp(α i+ 1 2 ) -1 D i+ 1 2 ∆ε i+ 1 2 f i .
And now, we can write

K i+ 1 2 = A i+ 1 2 f i+1 -B i+ 1 2 f i .
Summary To summarize we name scheme S 1 the discretized system (3.7) where the numerical fluxes are (for 1 ≤ i ≤ N -1) to be identically equal to zero. This strategy was carried by Langdon [START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF] and Decoster and Langdon [START_REF] Decoster | [END_REF] in the case of Coulombian interactions. We extend their method to the general case. We obtain (r.h.s.) = -1 2

K i+ 1 2 = E i+ 1 2 f i+ 1 2 + D i+ 1 2 (f i+1 -f i )/∆ε i+ 1 2 with E i+ 1 2 = - N -1 j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 , D i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 and f i+
N i=1 N j=1 ((Dφ) i+ 1 2 -(Dφ) j+ 1 2 ) g i+ 1 2 ,j+ 1 2 (f j+ 1 2 (Df ) i+ 1 2 -f i+ 1 2 (Df ) j+ 1 2 ) ∆ε i+ 1 2 ∆ε j+ 1 2 ,
(3.17) and

N i=1 c i ∂f i ∂t φ i = - N i=1 (Dφ) i+ 1 2 ∆ε i+ 1 2 N j=1 g i+ 1 2 ,j+ 1 2 (f j+ 1 2 (Df ) i+ 1 2 -f i+ 1 2 (Df ) j+ 1 2 ) ∆ε j+ 1 2 .
(3.18) Finally, we get the weak formulation of the semi-discretized model

df i dt = Q S 2 i 1 ≤ i ≤ N (3.19) with Q S 2 1 = K 3 2 /c 1 , Q S 2 i = ( K i+ 1 2 -K i- 1 2
)/c i for 2 ≤ i ≤ N , where

K i+ 1 2 = -f i+ 1 2 N j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 + (Df ) i+ 1 2 N j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 , 1 ≤ i ≤ N
(3.20) Like for the previous scheme we can write

E i+ 1 2 = - N j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 and D i+ 1 2 = N j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 , thus K i+ 1 2 = E i+ 1 2 f i+ 1 2 + D i+ 1 2 f i+1 -f i ∆ε i+ 1 2 .

Now we have

K N + 1 2 = E N + 1 2 f N + 1 2 + D N + 1 2 f N +1 -f N ∆ε N + 1 2 ,
as f N +1 = 0 the flux at the last point reads

K N + 1 2 = E N + 1 2 f N + 1 2 -D N + 1 2 f N ∆ε N + 1 2 . We enforce K N + 1 2 to vanish (to impose df N +1 dt = 0
). Thus we get the relation

E N + 1 2 f N + 1 2 -D N + 1 2 f N ∆ε N + 1 2 = 0 . (3.21)
Thanks to (3.21) the diffusion becomes

D i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + g i+ 1 2 ,N + 1 2 D N + 1 2 E N + 1 2 f N . (3.22)
Remark 7. We integrate by parts the drift term and get

E i+ 1 2 = N j=1 (g i+ 1 2 ,j+ 1 2 -g i+ 1 2 ,j-1 2 ) f j ∆ε j+ 1 2 .
We assume that {g i+ 1 2 ,j+ 1 2 } 1≤i≤N is an increasing sequence. Then, if the f i 's are positive,

E i+ 1 2 is positive too. In other hand, even if the f i+ 1 2 's are positive D i+ 1 2 can be negative if f N = 0. The discretization of f i+ 1 2
is still given by the Chang and Cooper average. Therefore we get

f i+ 1 2 = (1 -δ i+ 1 2 )f i+1 + δ i+ 1 2 f i , ∀ i ; 1 ≤ i ≤ N -1 , (3.23)
where

δ i+ 1 2 = h( α i+ 1 2 ) = 1 α i+ 1 2 - 1 (exp( α i+ 1 2 ) -1) , ∀ i ; 1 ≤ i ≤ N -1 , (3.24)
and with α i+

1 2 = E i+ 1 2 D i+ 1 2 ∆ε i+ 1 2
. Thanks to the boundary condition (3.21) we can extend the

definition to N + 1 2 by setting δ N + 1 2 = 1/ α N + 1 2
(see [START_REF] Decoster | [END_REF]). Now α i+ 1 2 reads

α i+ 1 2 =         - N j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 ((1 -δ j+ 1 2 )f j+1 + δ j+ 1 2 f j ) ∆ε j+ 1 2 + g i+ 1 2 ,N + 1 2 D N + 1 2 E N + 1 2 f N         ∆ε i+ 1 2 ,
As in the previous scheme we suppose that δ is the N-dimensional column vector with components

{ δ i+ 1 2 } 1≤i≤N . Thus we can write δ i+ 1 2 as a function of f and introduce { H i+ 1 2 } 1≤i≤N such that δ i+ 1 2 (f ) = H i+ 1 2 ( δ 3 2 , ..., δ k+ 1 2 , ..., δ N + 1 2 , f 1 , ..., f k , ... f N ) ,
in other words, to compute δ, the system to be solved is

δ(f ) = H( δ(f ), f ) , (3.25)
where f is a distribution function that has density ρ and energy ρE and H is a N-dimensional vector valued function. Here the problem rises from the fact that { D i+ 1 2 } 1≤i≤N can be undefined even if f is properly defined. We'll specify this more clearly in the case where potential are Coulombian or Maxwellian.

satisfies the relation

K i+ 1 2 = A i+ 1 2 f i+1 -B i+ 1 2 f i (3.26)
where

B i+ 1 2 = u( α i+ 1 2 ) D i+ 1 2 ∆ε i+ 1 2 , (3.27)
and

A i+ 1 2 = v( α i+ 1 2 ) D i+ 1 2 ∆ε i+ 1 2 , (3.28
)

with α i+ 1 2 = E i+ 1 2 D i+ 1 2 ∆ε i+ 1 2 , u(α) = α exp(α) -1 and v(α) = α exp(α) exp(α) -1 .
The proof is checked as in Proposition 1.

Summary. To summarize we identify as scheme S 2 the system (3.19) with boundary conditions (3.21) where the numerical fluxes are

K i+ 1 2 = E i+ 1 2 f i+ 1 2 + D i+ 1 2 (f i+1 -f i )/∆ε i+ 1 2 with E i+ 1 2 = - N j=1 g i+ 1 2 ,j+ 1 2 (Df ) j+ 1 2 ∆ε j+ 1 2
and

D i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + g i+ 1 2 ,N + 1 2 f N D N + 1 2 / E N + 1 2
. Here f i+ Remark 8. Comparing formally the drift and diffusion functional discretized respectively by S 1 and S 2 we get for

1 ≤ i ≤ N -1 E i+ 1 2 = E i+ 1 2 + f N g i+ 1 2 ,N + 1 2 and D i+ 1 2 = D i+ 1 2 + ∆ε N + 1 2 f N + 1 2 g i+ 1 2 ,N + 1 2 . (3.29)
Remark 9. Langdon [START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF] has previously developed a scheme for Coulombian interactions. This scheme is nothing else than S 2 reduced to the Coulombic potentials.

A particular case: the Maxwellian interactions

First we focus on the expression of the collision terms in the case where potentials are Maxwellian. The diffusion-drift ratio is constant and consequently the Chang and Cooper relation is reduced to a non-linear scalar equation. As a result we'll show that for the scheme S 2 the diffusion term can be undefined.

For Maxwellian interactions we have

g i+ 1 2 ,j+ 1 2 = ε 3 2 i+ 1 2 ε 3 2 j+ 1 2 . Then, E i+ 1 2
can be write by factorizing the term f i

E i+ 1 2 = ε 3 2 i+ 1 2 N -1 j=1 (ε 3 2 j+ 1 2 -ε 3 2 j-1 2 ) f j -f N ε 3 2 N -1 2 ε 3 2 i+ 1 2
and we get for the diffusion term

D i+ 1 2 = ε 3 2 i+ 1 2 N -1 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 .
Now we turn to the scheme S 2 and obtain in the Maxwellian case, assuming that f is a distribution function that has density ρ and energy ρE and that (3.21) holds,

               E i+ 1 2 = ε 3 2 i+ 1 2 N j=1 (ε 3 2 j+ 1 2 -ε 3 2 j-1 2 ) f j = 3 2 ε 3 2 i+ 1 2 N j=1 √ ε j ∆ε j f j = 3 2 ε i+ 1 2 ρ , D i+ 1 2 = ε 3 2 i+ 1 2 N j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 f j+ 1 2 = ε 3 2 i+ 1 2 N -1 j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 f j+ 1 2 + ε 3 2 i+ 1 2 ε 3 2 N + 1 2 f N D N + 1 2 E N + 1 2 , (3.30
) thus

E i+ 1 2 = E i+ 1 2 + ε 3 2 i+ 1 2 ε 3 2 N + 1 2 f N and D i+ 1 2 = D i+ 1 2 + ε 3 2 i+ 1 2 ε 3 2 N + 1 2 f N D N + 1 2 E N + 1 2 , (3.31)
as expected from (3.29).

Remark 10. We have (3.21) and (3.31) we can write

D i+ 1 2 /E i+ 1 2 = C for i; 1 ≤ i ≤ N -1 and D i+ 1 2 / E i+ 1 2 = C for i; 1 ≤ i ≤ N . Following
D i+ 1 2 E i+ 1 2 =         D i+ 1 2 E N + 1 2 D N + 1 2 + ε 3 2 i+ 1 2 ε 3 2 N + 1 2 f N E i+ 1 2 + ε 3 2 i+ 1 2 ε 3 2 N + 1 2 f N         D N + E N + 1 2 . As ( D i+ 1 2 / E i+ 1 2 ) = ( D N + 1 2 / E N + 1 2 ) we get C = C.
Now, we are interested by the expression of δ i+ 1 2 in the case where potentials are Maxwellian. Recall that with scheme S 1 we get

α i+ 1 2 = E i+ 1 2 /D i+ 1 2 ∆ε i+ 1 2 and h(α) = 1/α -1/(exp(α) -1).
Thus, thanks to (3.30) and (3.31) we can write

α i+ 1 2 = 3 2 ρ -f N ε 3 2 N + 1 2 N -1 j=1 ε 3 2 j+ 1 2 ∆ε j+ 1 2 f j+ 1 2 ∆ε i+ 1 2 . We set d(f ) = N -1 j=1 ε 3 2 j+ 1 2 ∆ε j+ 1 2 f j+ 1 2 therefore as f i+ 1 2
is given by the Chang and Cooper average

d(f ) = N -1 j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 ((1 -δ j+ 1 2 )f j+1 + δ j+ 1 2 f j ) ,
and introducing h(α i+ 1 2

) we get

d(f ) = N -1 j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 (h( 3 2 ρ -f N ε 3 2 N + 1 2 d(f ) ∆ε j+ 1 2 ) (f j -f j+1 ) + f j+1 ) . ( 3
N i=1 c i ∂f i ∂t φ i = - N -1 i=1 (Dφ) i+ 1 2 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 (f i+1 f j -f i f j+1 ),
which is entropic as shown in [START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF].

Following the same way we get for the scheme S 2

α i+ 1 2 = 3 2 ρ N -1 j=1 ε 3 2 j+ 1 2 ∆ε j+ 1 2 f j+ 1 2 + ε 3 2 N + 1 2 f N + 1 2 ∆ε N + 1 2 ∆ε i+ 1 2 .
We denote

d(f ) = N -1 j=1 ε 3 2 j+ 1 2 ∆ε j+ 1 2 f j+ 1 2 + ε 3 2 N + 1 2 f N + 1 2 ∆ε N + 1 2
where f i+ 1 2 is computed by the Chang and Cooper average. Thus

d(f ) = N -1 j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 ((1 -δ j+ 1 2 )f j+1 + δ j+ 1 2 f j ) + ε 3 2 N + 1 2 f N + 1 2 ∆ε N + 1 2
, so that, according that (3.21) holds and introducing h( α i+ 1

2

), we have

d(f ) = N -1 j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 (h( 3 2 ρ d(f ) ∆ε j+ 1 2 ) (f j -f j+1 ) + f j+1 ) + ε 3 2 N + 1 2 f N d(f ) 3 2 ρ . (3.33) Hence if 3 2 ρ = ε 3 2 N + 1 2 f N the term d(f ) is undefined. As the f i+ 1 2
's are positive we remark that d(f ) and

Y = ( 3 2 ρ -ε 3 2 N + 1 2 f N ) (3.34)
have the same sign. To find δ i+ 1 2

we have to solve the non-linear equation (3.33).

A particular case: the Coulombian interactions

In a first part we give the expression of drift and diffusion in the case where the potentials are Coulombian. For S 2 we find Langdon's scheme as planned. We point out that { D i+ 1 2 } 1≤i≤N (diffusion term computed by the scheme S 2 ) can be undefined. After that we carry on in offering an equivalent form of (3.13) (respectively (3.25)) to compute the interpolant of the Chang and Cooper average as in [START_REF] Decoster | [END_REF]. An advantage is that this form requires less operation to implement.

For Coulombian interactions we get

g i+ 1 2 ,j+ 1 2 = g(ε i+ 1 2 , ε j+ 1 2 ) = min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2
). This relation is inserted in (3.9) to obtain the drift term discretized by the scheme S 1

E i+ 1 2 = - N -1 j=1 ∆ε j+ 1 2 min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2 ) f j+1 -f j ∆ε j+ 1 2 for i ; 1 ≤ i ≤ N -1 ,
by factorizing the term f i we have

E i+ 1 2 = N -1 j=1 (min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2
) -min(ε

3 2 i+ 1 2 , ε 3 2 j-1 2
)) f j -f N min(ε 3 2 i+ 1 2 , ε 3 2 N - 1 2
) .

(3.35)

Now {ε i+ 1 2
} 1≤i≤N is an increasing sequence, therefore

E i+ 1 2 = i j=1 (ε 3 2 j+ 1 2 -ε 3 2 j-1 2 ) f j -f N min(ε 3 2 i+ 1 2 , ε 3 2 N -1 2 ) ∀ i ; 1 ≤ i ≤ N -1 . (3.36)
In the same way, we get for the diffusion term

D i+ 1 2 = N -1 j=1 ∆ε j+ 1 2 min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2 ) f j+ 1 2 ∀ i ; 1 ≤ i ≤ N -1 , (3.37)
or

             D i+ 1 2 = i j=1 ∆ε j+ 1 2 f j+ 1 2 ε 3 2 j+ 1 2 + ε 3 2 i+ 1 2 N -1 j=i+1 ∆ε j+ 1 2 f j+ 1 2 , 1 ≤ i ≤ N -2 D N -1 2 = N -1 j=1 ∆ε j+ 1 2 f j+ 1 2 ε 3 2 j+ 1 2 . ( 3
E i+ 1 2 = - N j=1 min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2 )(f j+1 -f j ) = N j=1
(min(ε

3 2 i+ 1 2 , ε 3 2 j+ 1 2
)-min(ε

3 2 i+ 1 2 , ε 3 2 j-1 2
))f j , (3.39) 20

E i+ 1 2 = i j=1 (ε 3 2 j+ 1 2 -ε 3 2 j-1 2 ) f j and E N -1 2 = 3 2 ρ -(ε 3 2 N + 1 2 -ε 3 2 N -1 2 ) f N , (3.40)
and

D i+ 1 2 = N -1 j=1 ∆ε j+ 1 2 min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2 )f j+ 1 2 + ∆ε N + 1 2 f N + 1 2 ε 3 2 i+ 1 2 , (3.41)
or,

             D i+ 1 2 = i j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ε 3 2 i+ 1 2 N -1 j=i+1 f j+ 1 2 ∆ε j+ 1 2 + ∆ε N + 1 2 f N + 1 2 ε 3 2 i+ 1 2 , 1 ≤ i ≤ N -2 D i+ 1 2 = N -1 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ∆ε N + 1 2 f N + 1 2 ε 3 2 i+ 1 2 , i = N -1 or i = N .
Note that contrary to scheme S 1 the drift and diffusion terms computed by scheme S 2 are not consistent approximations of E(f ) and D(f ) when their integration domain is reduced to a bounded domain [0, E] in the variable ε (see equations (2.13) and (2.14)). To compare to the Langdon scheme, we prefer to write

E i+ 1 2 , D i+ 1 2 , E i+ 1 2 and D i+ 1 2
with the velocity variable when it is necessary. Thanks to (3.2) we have

E i+ 1 2 = 3 i j=1 (v 2 ∆v) j f j -f N min(ε 3 2 i+ 1 2 , ε 3 2 N -1 2
) .

(3.42)

Seeing ε 1 2 = 0 we can write ε

3 2 i+ 1 2 = i j=1 (ε 3 2 j+ 1 2 -ε 3 2 j- 1 2
) and from (3.2) we get ε

3 2 i+ 1 2 = 3 i j=1 (v 2 ∆v) j . So              D i+ 1 2 = 3 i j=1 ∆ε j+ 1 2 f j+ 1 2 j k=1 (v 2 ∆v) k + 3 N -1 j=i+1 ∆ε j+ 1 2 f j+ 1 2 i k=1 (v 2 ∆v) k , 1 ≤ i ≤ N -2 D N -1 2 = 3 N -1 j=1 ∆ε j+ 1 2 f j+ 1 2 j k=1 (v 2 ∆v) k
By factorizing the term (v 2 ∆v) k we write

D i+ 1 2 = 3 i j=1 (v 2 ∆v) j N -1 k=j f k+ 1 2 ∆ε k+ 1 2 1 ≤ i ≤ N -1 . (3.43)
For the scheme S 2 we obviously find

E i+ 1 2 = 3 i j=1 (v 2 ∆v) j f j and D i+ 1 2 = 3 i j=1 (v 2 ∆v) j N k=j f k+ 1 2 ∆ε k+ 1 2
as in Langdon's scheme [START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF]. That is, following (3.42) and (3.43), for 1

≤ i ≤ N -1 E i+ 1 2 = E i+ 1 2 -f N min(ε 3 2 i+ 1 2 , ε 3 2 N -1 2 ) and D i+ 1 2 = D i+ 1 2 -f N + 1 2 ε 3 2 i+ 1 2 ∆ε N + 1 2 . ( 3
E N + 1 2 = E N -1 2 + ε 3 2 N + 1 2 f N and D N + 1 2 = D N -1 2 + ε 3 2 N + 1 2 ∆ε N + 1 2 f N + 1 2 = D N -1 2 + ε 3 2 N + 1 2 f N D N + E N + , consequently D N + 1 2 E N + 1 2 = D N + 1 2 E N + 1 2 D N -1 2 E N + 1 2 D N + 1 2 + ε 3 2 N + 1 2 f N E N -1 2 + ε 3 2 N + 1 2 f N .
This leads to the relation

D N + 1 2 E N + 1 2 = D N -1 2 E N -1 2 .
Using (3.21) and (3.44) again, we write

E N -1 2 = E N -1 2 + ε 3 2 N -1 2 f N and D N -1 2 = D N -1 2 + ε 3 2 N -1 2 f N D N + 1 2 E N + 1 2 , hence D N + 1 2 E N + 1 2 = D N -1 2 E N -1 2 .
(3.45)

Decoster [START_REF] Decoster | [END_REF] shows that the relation (3.45) is sufficient to conserve energy and imposes it as boundary condition. Finally this is not necessary in view that (3.45) is held in the discretization.

Here the diffusion term discretized by S 2 reads, thanks to (3.21) and (3.45) 's are positive. The positivity condition f N ε

D i+ 1 2 = N -1 j=1 ∆ε j+ 1 2 min(ε 3 2 i+ 1 2 , ε 3 2 j+ 1 2 )f j+ 1 2 + ε 3 2 i+ 1 2 f N D N -1 2 E N -1 2 , (3.46) therefore D N -1 2 = N -1 j=1 ∆ε j+ 1 2 ε 3 2 j+ 1 2 f j+ 1 2 + ε 3 2 N -1 2 f N D N -1 2 E N -1 2 . (3.47) Now if f N ε 3 2 N -1 2 = E N -1 2 the term D N -
Y = ( E N -1 2 -f N ε 3 2 N -1 2 ) = 3 2 ρ -f N ε 3 2 N + 1 2 (3.
3 2 N -1 2 < E N -1 2
is the same condition given by Decoster to avoid what he calls "abnormal discretization" [START_REF] Decoster | [END_REF]. Note that if D N -1 2 is undefined the system (3.25) doesn't have any solution.

In this particular "Coulombian" case an other way to compute the Chang and Cooper average interpolant δ i+ 1 2 is to consider equation (2.18) as in [START_REF] Decoster | [END_REF] in place of system (3.13). Thus, we have to solve an elliptic equation, namely

- 2 3 ∂ ∂ε 1 √ ε ∂ ∂ε (D(f, ε)) = f (ε) in R + , with boundary condition D(f, 0) = 0. The solution is D(f ) = ε 0 ε 3 2 f (ε )dε + ∞ ε f (ε )dε .
The non-linear system (3.13) which solution is δ(f ) is then replaced by a non-linear system on D(f ).

Proposition 3. The diffusion term discretized by S 1 is solution of the following system

M D = f δ (3.49)
where D (respectively f δ ) is the (N -1)-dimensional column vector with components

{D i+ 1 2 } 1≤i≤N -1 (respectively {∆ε i+ 1 2 f i+ 1 2 } 1≤i≤N -1
) and M is the (N -1)×(N -1) tridiagonal matrix with components

M 1,2 = -(ε 3 2 5 2 -ε 3 2 3 2 ) -1 , M 1,1 = -M 1,2 + (ε 3 2 3 2 ) -1 M i,i-1 = -(ε 3 2 i+ 3 2 -ε 3 2 i+ 1 2 ) -1 , M i,i+1 = -(ε 3 2 i+ 1 2 -ε 3 2 i-1 2 ) -1 , M i,i = -(M i,i-1 + M i,i+1 ) , M N -1 2 ,N -3 2 = -(ε 3 2 N -1 2 -ε 3 2 N -3 2 ) -1 , M N -1 2 ,N -1 2 = -M N -1 2 ,N -3 2 , Proof. From (3.38) we get for i; 1 ≤ i ≤ N -2 D i+ 1 2 = i j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ε 3 2 i+ 1 2 N -1 j=i+1 f j+ 1 2 ∆ε j+ 1 2 ,
and for i;

2 ≤ i ≤ N -1 D i-1 2 = i-1 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ε 3 2 i-1 2 N -1 j=i f j+ 1 2 ∆ε j+ 1 2 .
By subtracting the second equation from the first one we get

D i+ 1 2 -D i-1 2 = ε 3 2 i+ 1 2 f i+ 1 2 ∆ε i+ 1 2 + (ε 3 2 i+ 1 2 -ε 3 2 i-1 2 ) N -1 j=i+1 f j+ 1 2 ∆ε j+ 1 2 -ε 3 2 i-1 2 f i+ ∆ε i+ 1 2 . Thus D i+ 1 2 -D i-1 2 ε 3 2 i+ 1 2 -ε 3 2 i-1 2 = f i+ 1 2 ∆ε i+ 1 2 + N -1 j=i+1 f j+ 1 2 ∆ε j+ 1 2 .
In the same way we have the relation

D i+ 3 2 -D i+ 1 2 ε 3 2 i+ 3 2 -ε 3 2 i+ 1 2 = f i+ 3 2 ∆ε i+ 3 2 + N -1 j=i+2 f j+ 1 2 ∆ε j+ 1 2 , hence D i+ 3 2 -D i+ 1 2 ε 3 2 i+ 3 2 -ε 3 2 i+ 1 2 - D i+ 1 2 -D i-1 2 ε 3 2 i+ 1 2 -ε 3 2 i-1 2 = -f i+ 1 2 ∆ε i+ 1 2 , ∀i ; 2 ≤ i ≤ N -2 . (3.50)
Now we have to produce the relationship associated to i = N -1 and i = 2. First, from (3.37) we get relations

D N -1 2 = N -1 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 and D N -3 2 = N -2 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ε 3 2 N -3 2 f N - ∆ε N -1 2 ,
and so on

D N -1 2 -D N -3 2 ε 3 2 N -1 2 -ε 3 2 N -3 2 = f N -1 2 ∆ε N -1 2 . (3.51) Furthermore D 3 2 = ε 3 2 3 2 f 3 2 ∆ε 3 2 + ε 3 2 3 2 N -1 j=2 f j+ 1 2 ∆ε j+ 1 2 .
We sum the (N -3) relations (3.50) from i = 2 to i = N -2 and (3.51). Therefore

D 3 2 = ε 3 2 3 2 f 3 2 ∆ε 3 2 + ε 3 2 3 2 D 5 2 -D 3 2 ε 3 2 5 2 -ε 3 2 3 2
.

This completes the proof.

Remark 12. The system (3.49) can be solved by a Newton method in only O(N ) operations at each step.

We can make also the following remark

Remark 13. M is a M-matrix. Actually, according that {ε i+ 1 2 } 1≤i≤N -1 is an increasing sequence M i,i-1 and M i,i+1 are negative. Furthermore

M i,i = -(M i,i-1 + M i,i+1 ) then M i,i > 0. As M is non-singular M is a M-matrix (that means (M -1 v, v) > 0; ∀v > 0). So, if {f i+ 1 2 } 1≤i≤N -1 is positive the D i+ 1 2
's remain positive. This confirms Remark 6. We'll see latter that this property is unavoidable to guarantee the positivity of the scheme S 1 .

We turn now to the scheme S 2 (always in the particular case where potential are Coulombian) and produce the following proposition. Proposition 4. Let f a distribution function that has density ρ and energy ρE. We assume

that f N ε 3 2 N -1 2 / E N -1 2 = 1.
The diffusion term discretized by S 2 is solution of the following system

M D = f e δ where D (respectively f e δ ) is the (N -1)-dimensional column vector with components { D i+ 1 2 } 1≤i≤N -1 (respectively {∆ε i+ 1 2 f i+ 1 2 } 1≤i≤N -1
) and M is the (N -1)×(N -1) tridiagonal matrix with components

M 1,2 = -(ε 3 2 5 2 -ε 3 2 3 2 ) -1 , M 1,1 = -M 1,2 + (ε 3 2 3 2 ) -1 M i,i-1 = -(ε 3 2 i+ 3 2 -ε 3 2 i+ 1 2 ) -1 , M i,i+1 = -(ε 3 2 i+ 1 2 -ε 3 2 i-1 2 ) -1 , M i,i = -( M i,i-1 + M i,i+1 ) , M N -1 2 ,N -3 2 = -(ε 3 2 N -1 2 -ε 3 2 N -3 2 ) -1 , M N -1 2 ,N -1 2 = -M N -1 2 ,N -3 2 - f N E N -1 2 .
Proof. As in the previous proof we easily check that

D i+ 3 2 -D i+ 1 2 ε 3 2 i+ 3 2 -ε 3 2 i+ 1 2 - D i+ 1 2 -D i-1 2 ε 3 2 i+ 1 2 -ε 3 2 i-1 2 = -f i+ 1 2 ∆ε i+ 1 2 , 2 ≤ i ≤ N -2 . (3.52)
We look what's happening at i = N -1. From the discretization of the diffusion term we get

D N -1 2 = N -1 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ε 3 2 N -1 2 f N + 1 2 ∆ε N + 1 2
, and

D N -3 2 = N -2 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 + ε 3 2 N -3 2 f N -1 2 ∆ε N -1 2 + ε 3 2 N -3 2 f N + 1 2 ∆ε N + 1 2 . Therefore D N -1 2 -D N -3 2 ε 3 2 N -1 2 -ε 3 2 N -3 2 = f N -1 2 ∆ε N -1 2 + f N + 1 2 ∆ε N + 1 2 .
But following (3.21) and (3.45) we can write

f N + 1 2 ∆ε N + 1 2 = f N D N -1 2 / E N -1 2
, thus Note that M is a M-matrix if and only if

D N -1 2 ( 1 ε 3 2 N -1 2 -ε 3 2 N -3 2 - f N E N -1 2 ) - D N -3 2 ε 3 2 N -1 2 -ε 3 2 N -3 2 = f N -1 2 ∆ε N -1 2 . ( 3
f N < E N -1 2 (ε 3 2 N -1 2 -ε 3 2 N -3 2 ) -1 and f N ε 3 2 N -1 2 / E N -1 2 = 1.

Mass and energy conservations

Proposition 5. The schemes S 1 and S 2 are conservative in density and energy.

These conservation properties can be easily checked by taking {φ

i } 1≤i≤N = {(1, ε i ) t } 1≤i≤N in (3.5
) and (3.17) respectively.

Equilibrium solution

Recall that M = exp(-βε) is the Maxwellian equilibrium solution of the equation (2.1) that has the same density and energy as f (ε, 0) = f 0 (ε). We note {M i } 1≤i≤N = {exp(-βε i )} 1≤i≤N the approximation of {M (ε i )} 1≤i≤N on the energy grid. In this section, we focus on the thermodynamical equilibrium state. We show in particular that the schemes S 1 and S 2 preserve the equilibrium state when is reached; it means that the numerical fluxes are equal to zero when the distribution f is equal to the equilibrium solution. First we consider the scheme S 1 .

Proposition 6. We assume that M i+ 1 2 is computed by Chang and Cooper average formula (3.11) where

δ i+ 1 2 = 1 β∆ε i+ 1 2 - 1 exp(β∆ε i+ 1 2 ) -1 ∀i ; 1 ≤ i ≤ N -1 . (3.54)
Thus we have

M i+ 1 2 = - M i+1 -M i β∆ε i+ 1 2 and E i+ 1 2 D i+ 1 2 = β for all i such that 1 ≤ i ≤ N -1 where E i+ 1 2
and D i+ 1 2 are computed by the scheme S 1 (see (3.9)).

Proof. First we remark that M i+1 = exp(-β∆ε i+ 1 2

)M i thus we have

M i+1 -M i = -M i+1 (exp(-β∆ε i+ 1 2 ) -1) ,
Now from (3.11) and (3.54) we get

M i+ 1 2 = 1 β∆ε i+ 1 2 - 1 exp(-β∆ε i+ 1 2 ) -1 (M i -M i+1 ) + M i+1 . Therefore M i+ 1 2 = ( 1 β∆ε i+ 1 2 + M i+1 M i+1 -M i ) (M i -M i+1 ) + M i+1 , and
M i+ 1 2 = - M i+1 -M i β∆ε i+ 1 2 . (3.55)
Thanks to (3.55) the drift term reads at equilibrium

E i+ 1 2 = - N -1 j=1 g i+ 1 2 ,j+ 1 2 (M j+1 -M j ) = N -1 j=1 g i+ 1 2 ,j+ 1 2 M j+ 1 2 β∆ε j+ 1 2 = βD i+ 1 2 .
That completes the proof.

Note that the Proposition 6 amounts to say that problem (3.13) has a solution at equilibrium given by relation (3.54). We are now going to show that at equilibrium this solution is unique. By the Brouwer fixed point theorem we saw that H has a fixed point in [0, 1] N -1 .

To prove the uniqueness of solution we just have to show that the mapping H is contracting, i.e., H(x) -H(y) ≤ L x -y for some norm, Lipschitz constant L < 1, and all x, y in [0, 1] N -1 . Indeed, if H has two fixed points δ and δ we get H(δ) -H(δ ) < δ -δ and that is conflicting. We introduce the matrix a(δ) = (a i,j ) = ( ∂H i+ 1 2 (δ)

∂δ j+ 1 2
) and the matrix norm induced by the maximum norm a = max i N -1 j=1 | a i,j |. At equilibrium we have

a i,j = -h (α i+ 1 2 ) E i+ 1 2 D 2 i+ 1 2 ∆ε i+ 1 2 g i+ 1 2 ,j+ 1 2 (M j -M j+1 )∆ε j+ 1 2 .
Note that as M is a Maxwellian, M is monotone thus M j -M j+1 ≥ 0 and

E i+ 1 2 ≥ 0. Therefore N -1 j=1 | a i,j |=| h (α i+ 1 2 ) | E i+ 1 2 D 2 i+ 1 2 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 (M j -M j+1 )∆ε j+ 1 2 .
We denote ∆ε = max i (∆ε i+ 1 2

) and we get

N -1 j=1 | a i,j | ≤ | h (α i+ 1 2 ) | α 2 i+ 1 2 ∆ε ∆ε i+ 1 2 . Now let g(x) = h (x)x 2
. By calculating g it's easy to see that the function g is decaying over [0, +∞[ and g(0) = 0, g(+∞) = -1. Therefore, since x is bounded we have

| g(x) |< 1. By definition α i+ 1 2 reads α i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 M j+1 (exp(β∆ε j+ 1 2 ) -1) N -1 j=1 g i+ 1 2 ,j+ 1 2 ((1 -δ j+ 1 2 )M j+1 + δ j+ 1 2 M j ) ∆ε i+ 1 2 . Because E i+ 1 2 and D i+ 1 2 are positive, δ i+ 1 2 lies in [0, 1/2]. Hence α i+ 1 2 ≤ 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 M j+1 (exp(β∆ε j+ 1 2 ) -1) N -1 j=1 g i+ 1 2 ,j+ 1 2 M j+1 ∆ε j+ 1 2 ∆ε i+ 1 2 .
But (exp(x) -1)/x is an increasing function, consequently we have We turn now to scheme the S 2 .

α i+ 1 2 ≤ 2 exp(β∆ε) -1 ∆ε ∆ε i+ 1 2 , and | g(α i+ 1 2 ) |< 1. It is clear that if {∆ε i+ 1 2 } 1≤i≤N is a uniform sequence then N -1 j=1 | a i,
Proposition 7. We assume that M i+ 1 2 is computed by the Chang and Cooper average formula (3.23) where

δ i+ 1 2 = 1 β∆ε i+ 1 2 - 1 exp(β∆ε i+ 1 2 ) -1 ∀i ; 1 ≤ i ≤ N -1 .

Thus we have

E i+ 1 2 D i+ 1 2 = β for all i such that 1 ≤ i ≤ N where E i+ 1 2
and D i+ 1 2 are calculated with the scheme S 2 .

Proof. According to (3.29) and thanks to (3.21) we can write

D N + 1 2 E N + 1 2 = D N + 1 2 E N + 1 2 D N + 1 2 E N + 1 2 D N + 1 2 + g N + 1 2 ,N + 1 2 f N E N + 1 2 + g N + 1 2 ,N + 1 2 f N . Therefore D N + 1 2 E N + 1 2 = D N + 1 2 E N + 1 2 = 1 β . (3.56)
Once again we get

D i+ 1 2 E i+ 1 2 = D N + 1 2 E N + 1 2 D i+ 1 2 E N + 1 2 D N + 1 2 + g i+ 1 2 ,N + 1 2 f N E i+ 1 2 + g i+ 1 2 ,N + 1 2 f N ,
and consequently

D i+ 1 2 E i+ 1 2 = D i+ 1 2 E i+ 1 2 = 1 β
. This ends the proof.

Positivity

If equation (2.1) is linear (in the sense that E(f ) and D(f ) are independent of f ), Chang and Cooper [START_REF] Chang | A practical difference scheme for Fokker-Planck equations[END_REF] have shown that their scheme is positive; in other words, if initial value f 0 (ε) is non-negative then the value of f (ε, t) at each successive time step remains non-negative. Afterwards, Larsen et al. [START_REF] Larsen | Discretization methods for one-dimensional Fokker-Planck operators[END_REF] have pointed out that if (2.1) is non-linear (in particular if there exists non-linearity in energy variable) the Chang and Cooper method generally loses the property that numerical solutions are guaranteed to be positive. We are interested in the case where E(f ) and D(f ) depend on the distribution function. We show that S 1 is a positive scheme while S 2 can produce negative solution. Now we focus on the scheme S 1 .

Proposition 8. We assume that for all distribution function f = {f i } 1≤i≤N there exists

δ(f ) = {δ i+ 1 2 (f )} 1≤i≤N -1
that is C 1 and given by (3.12). Then (3.7) has a solution. Furthermore, if f i (t = 0) is non-negative, then f i (t) is non-negative for t > 0.

Proof. In the first place as δ(f ) is C 1 we are able to use Cauchy-Lipschitz's theorem. So that (3.7) has a unique solution. From the Proposition 1 we get

K i+ 1 2 = A i+ 1 2 f i+1 -B i+ 1 2 f i , with A i+ 1 2 = v(α i+ 1 2 ) D i+ 1 2 ∆ε i+ 1 2 and B i+ 1 2 = u(α i+ 1 2 ) D i+ 1 2 ∆ε i+ 1 2
and where

u(α) = α exp(α) -1 and v(α) = α exp(α) exp(α) -1 .
The function u(x) is positive, decreasing from -∞ to 0, while v(x) is positive too, increasing from 0 to +∞. To prove the positivity we use Lemma 1. Thus all we need to do is showing that the A i+ 1 are bounded. At first we assume that f is a distribution function that has density ρ and energy ρE thus we get

f j ≤ ρ max i c i , ∀j; 1 ≤ j ≤ N -1.
So that we can bound the drift term as following

| E i+ 1 2 |=| N -1 j=1 g i+ 1 2 ,j+ 1 2 (f j -f j+1 ) |≤ ρ max i c i 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 ≤ ρ max i c i 2(N -1)g N -1 2 ,N -1 2 .
Consequently, there exists

C > 0 such that | E i+ 1 2 |< C.
In the same way we get an upperbound for the diffusion coefficients

D i+ 1 2 ≤ ρ max i c i 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 ∆ε j+ 1 2 .
To continue we distinguish two cases. First we suppose that | α i+ 1 2

|≥ 1. Then we get is also bounded. Thus, thanks to Lemma 1 this ends the proof.

A i+ 1 2 = exp(α i+ 1 2 ) exp(α i+ 1 2 ) -1 E i+ 1 2 and B i+ 1 2 = 1 exp(α i+ 1 2 ) -1 E i+ 1 2 . And even if | α i+ 1 2 |→ ∞, A i+
Concerning the scheme S 2 we have shown that the discretized diffusion term can become negative. This loss of positivity prevent us to provide an efficient result for the positivity of the scheme considering that our proof is based on the D i+ 1 2 's positivity. All we can do is supplying an example of contrary in the case where potentials are Maxwellian. Let f a distribution function that as density ρ = N i=1 c i f i and energy ρE = N i=1 c i ε i f i . We assume that f = aδ ε 1 + bδ ε N where a and b are constants to define and the δ ε i 's are Dirac functions.

We assure also that the integration domain is sufficient to check that (

3 2 T = ρE ρ < ε N = E). Now we have ρ = c 1 a + c N b and ρE = c N ε N b (because ε 1 = 0). Thus b = ρE/(c N ε N ) and a = (ρ -ρE/ε N )/c 1 .
Note that a and b are non-negative. From (3.19) we have

c N -1 df N -1 dt = K N -1 2 -K N -3 2 = A N -1 2 f N -B N -1 2 f N -1 -( A N -3 2 f N -1 -B N -3 2 f N -2 ) , therefore df N -1 dt = A N -1 2 b/c N -1 . Now thanks to (3.28), A N -1 2 and D N -1 2
have the same sign.

Taking again what we saw previously (equation (3.33)) we write

D N -1 2 ε 3 2 N -1 2 (1 -ε 3 2 N + 1 2 f N 3 2 ρ ) = N -1 j=1 ε 3 2 j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 . Subsequently 1 -ε 3 2 N + 1 2 f N 3 2 ρ = 1 - ρE 3 2 ρ ε 3 2 N + 1 2 c N ε N = 1 -T ε 3 2 N + 1 2 ε N 1 √ ε N ∆ε N . If ∆ε N -→ 0 it is clear that D N -1 2
becomes negative. Therefore df N -1 dt < 0 and consequently f N -1 gets negative too. For the Coulombian case we can take the same example of contrary to show the non-positivity of the scheme S 2 .

We can be more precise with S 2 in the case where potentials are Maxwellian. Let us recall that

we denoted Y = 3 2 ρ -ε 3 2 N + 1 2 f N (see formula (3.34)).
As usual, f is a non-negative distribution that has non-zero density ρ and energy ρE such that T = 2/3 ρE ρ < T max . These functions generate a compact subset of

R N : D = {f ∈ R N ; ∀i, f i ≥ 0, i c i f i = ρ and i c i ε i f i = ρE}.
≤ δ i+ 1 2 ≤ 1 2 if Y > 0 and 1 2 ≤ δ i+ 1 2 ≤ 1 if Y < 0.
To begin, we assume that Y > 0. Thus, thanks to (3.33) we get the inequality

Y D N -1 2 /( 3 2 ρ) = ε 3 2 N -1 2 d(f )(1- ε 3 2 N + 1 2 f N 3 2 ρ ) = N -1 j=1 g N -1 2 ,j+ 1 2 f j+ 1 2 ∆ε j+ 1 2 ≥ N -1 j=1 g N -1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2 . The functional N -1 j=1 g N -1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2
is linear continuous on D thus achieves its minimum on this compact subset. This minimum is strictly positive because if it is equal to zero, that implies that all the mass will be concentrated at ε 1 , consequently the energy would be zero and we have supposed that the energy is strictly positive. On the other hand, if Y < 0, (3.33) gives

Y D N -1 2 /( 3 2 ρ) ≥ N -1 j=1 g N -1 2 ,j+ 1 2 f j ∆ε j+ 1 2 ,
and the right-hand side admits a strictly positive minimum (if not, all the mass will be concentrated at ε N and we have supposed that we only consider distribution function such that 0 ≤ T ≤ T max ). This shows that for the diffusion coefficient

D N -1 2
we ge lim

Y →±0 D N -1 2
→ ±∞. Since all the diffusion coefficient are proportionals and have the same sign, one has for all i: lim

Y →±0 D i+ 1 2 → ±∞.
In other words D i+ 1 2 behaves like 1 Y . This result says also that the scheme S 2 is not consistent with the truncated FPL equation defined by (3.1) but only with (2.3).

In the Coulombian case, the situation is more complicated. But we are able to construct a sequence of positive distribution functions f ν , with constant mass and energy, which have f 0 as limit as ν → ±0 and such that:

• f 0 verifies Y = 0, where Y is given by (3.48),

• there exist an index i 0 and a constant a such that min(f i 0 , f i 0 +1 ) ≥ a > 0,

• f 0 N -1 ≥ bf 0 N , with b > 2.
We begin by the construction of f 0 . First we assume that f 0 lies in D and we set

1 b f 0 N -1 = f 0 N = ε 3 2 N + 1 2 / 3 2 ρ.
The mass and the energy due to these two points are equal respectively to ∆ρ = ε

3 2 N + 1 2 (bc N -1 + c N )/ 3 2 ρ and ∆ρE = ε 3 2 N + 1 2 (bε N -1 c N -1 + ε N c N )/ 3 2 ρ.
Since c i ∆ε i √ ε i , for ∆ε N and ∆ε N -1 sufficiently small we have ρ-∆ρ > 0 and ρE-∆ρE > 0. Then we choose on the N -2 first points a positive distribution with mass ρ -∆ρ and with energy ρE -∆ρE, by example a Maxwellian having this mass and this energy. The resulting distribution f 0 verifies all the three points above.

The functional Y = 3 2 ρ -ε

3 2 N + 1 2 f N is continuous on the compact set D and {0} ∈ o Im(Y (D)),
since for the maxwellian that has density ρ and energy ρE, Y > 0 and for the example of contrary above Y < 0. Thus there exist sequences f ν ∈ D such that lim ν→0 f ν = f 0 with Y ν ≥ 0 or Y ν ≤ 0 and with Y ν → 0. For ν sufficiently small we have, by example, min(f ν i 0 , f ν i 0 +1 ) ≥ a/2 and f 0 N -1 > b/2bf 0 N . For this sequence Y ν > 0 we use the definition of the diffusion coefficients

Y ν D ν N -1 2 ≥ a 2 j g N -1 2 ,j+ 1 2 ∆ε j+ 1 2 3 2 ρ. Thus if Y ν ≥ 0 , lim ν→0 D ν N -1 2 = +∞ or if Y ν ≤ 0 , lim ν→0 D ν N -1 2 = -∞.
Since the definition of Y is the same in the Coulombian and in the Maxwellian cases, the above example works also in the Maxwellian case. We have seen that in the domain Y > 0 the diffusion coefficients are positive. And naturally the question that arise is: is this domain stable by (3.19)? it seems that this is not true: for the sequence f ν constructed above one can verify easily using (3.26) that df ν N dt tends to +∞ as Y ν tends to +0 that is dY ν dt tends to -∞ as Y ν tends to +0. This suggests that one can attain the boundary Y = 0 in finite time. And for Y = 0, as we have seen, the scheme S 2 is not defined.

Despite its deficiencies for the positivity one could use the scheme S 2 : a solution to avoid the lack of positivity would be to take the domain of computation sufficiently large in order to have the tail of the distribution function sufficiently small and thus Y We begin by given some remarks about the above analysis.

Remark 14. It is easy to check that the above results about positivity, energy conservation and equilibrium states for S 1 and S 2 don't depend on the definition of the coefficients c i , see

, or on our choice for the value of ε 1 .

Remark 15. For S 1 (and S 2 if the scheme is positive) if one takes the Chang and Cooper coefficients δ i+ 1 2 at equilibrium, that is given by (3.54), or freeze the value given by the initial data, it is not difficult to see that if it is done in the expression of the fluxes given by (3.10) the scheme is still positive but no more conservative in energy and if it is done in (3.14) the scheme is still conservative in energy but no more positive. The algebra that permits to pass from the form (3.10) to the form (3.14) is valid if and only if the coefficients verify at all times the relation (3.13). This remark is still valid if one does not calculate exactly the Chang and Cooper coefficients by an iterative method. Remark 16. The scheme S 2 is in fact not consistent with the truncated FPL equation (3.1), but only with (2.1).

Remark 17. We don't have any H-theorem for S 1 (and S 2 ) except for S 1 in the case where potentials are maxwellian and the energy grid is uniform, since in this case the scheme reduces to the entropic scheme of Berezin-Pekker. Remark 18. The scheme S 2 is not positive but it would not be a surprise if implicit time discretization leads to positive solutions. Actually it is well known that implicit time discretization adds some kind of numerical "viscosity" which "stabilizes" algorithms.

We are also interesting in the Epperlein velocity space discretization for solving the FPLE for Coulombian collisions [START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF] since it is based on the method of Chang and Cooper. This method is used in the code SPARK developed by Epperlein [START_REF] Epperlein | Fokker-Planck modelling of electrons transport in laser-produced plasmas[END_REF] or by Kingham and Bell in the recent code IMPACT [START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF]. As in the present paper the collisions coefficients are defined by Langdon. To difference the equation in velocity space, the author applies the Chang and Cooper approach. In his paper Epperlein develops a fully implicit finite-difference method for solving FPL equation and claims that his approach conserves both energy and number density exactly. Without talking about time discretization we just focus on his velocity space discretization and properties. In particular the Epperlein scheme differs on our scheme only by boundary conditions at ε = E. Actually, the author assumes that the diffusion term and the distribution function vanish at ε N + 1 2 and that the distribution function vanishes at ε N +1 . Thus the numerical flux is zero and consequently the density is conserved. In fact the Epperlein velocity discretization is nothing else than a miscellany of what we called scheme S 1 and scheme S 2 : the drift term is the one discretized by scheme S 2 , formula (3.39), and the diffusion term the one computed by S 1 , formula (3.37). Now, we can prove that the Epperlein velocity space discretization produces non-negative solution in consideration that the diffusion term is non-negative. Whereas we raise a doubt concerning the energy conservation. Effectively, assuming that diffusion vanishes at boundary while drift not, breaks the symmetry in the right-hand side of the semi-discretized equation. And the energy conservation proof leans on the symmetry between ε and ε (or v and v ). At last (3.56) fails and this relation is indispensable to check that the Chang and Cooper type scheme preserves the equilibrium state when is reached. As for S 2 this scheme is not consistent with the truncated FPL equation (3.1), but only with (2.1).

Remark 19. We can remark that these three Chang and Cooper type schemes are equivalent when the integrating domain is not bounded. The only difference between these schemes lies in the manner to treat boundary conditions when the integrating domain is reduced to [0, E].

Alternative schemes

The Chang and Cooper method is not the only way to provide positive, conservative and equilibrium states preserving schemes for the FPL equation. As we have seen above it is a very complicated scheme for this equation, not always positive or conservative in energy, depending on the boundary condition taken at the end of the domain of computation. The schemes we'll propose in the next sections share also these properties and are simpler.

Equilibrium scheme (scheme S 3 )

The first scheme is based on the work of Larsen et al. [START_REF] Larsen | Discretization methods for one-dimensional Fokker-Planck operators[END_REF]. Their work is for linear and nonlinear (in the sense that collisions terms are non-linear in energy variables) Fokker-Planck equations. But they do not consider non-linearity as in the Landau equation, that is drift and diffusion coefficients are functionnals of the distribution functions. One of the two main ideas exposed in their paper to preserve Maxwellian states is to remark that ε-derivative that appears in the flux can be rewritten as

∂ ∂ε = -βy ∂ ∂y
where y = exp(-βε). Thus it becomes easy to preserve Maxwellians in the weak form of the Fokker-Planck equation. Following these authors, in the case of the FPL equation, we set

∆ε i+ 1 2 = - (∆M ) i+ 1 2 βM i+1 = exp(β∆ε i+ 1 2 ) -1 β ∆ε i+ 1 2 (3.57)
Thus we consider the following approximation of the weak symmetrized form of the problem

N i=1 c i ∂f i ∂t φ i = - 1 2 N -1 i=1 N -1 j=1 ∆φ i+ 1 2 ∆ε i+ 1 2 - ∆φ j+ 1 2 ∆ε j+ 1 2 g i+ 1 2 ,j+ 1 2 f j+ 1 2 ∆f i+ 1 2 ∆ε i+ 1 2 -f i+ 1 2 ∆f j+ 1 2 ∆ε j+ 1 2 ∆ε j+ 1 2 ∆ε i+ 1 2 , ( 3
.58) with the approximations f i+ 1 2 taken decentered and defined by f i+ 1 2 = f i+1 . As in the previous section we obtain the system of ordinary equation for the approximation of the FPL equation

df i dt = Q S 3 i 1 ≤ i ≤ N (3.59)
where

Q S 3 1 = K 3 2 /c 1 , Q S 3 i = (K i+ 1 2 -K i-1 2 )/c i for 2 ≤ i ≤ N -1 and Q S 3 N = -K N -1 2 /c N with the numerical flux K i+ 1 2 = ∆ε i+ 1 2 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 (f j+1 ∆f i+ 1 2 ∆ε i+ 1 2 -f i+1 ∆f j+ 1 2 ∆ε j+ 1 2 )∆ε j+ 1 2 for i; 1 ≤ i ≤ N -1
. By factorizing the terms f i in the last sum, the numerical flux reads

K i+ 1 2 = ∆ε i+ 1 2 ∆ε i+ 1 2 ( ∆f i+ 1 2 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2 + f i+1 N -1 j=1 (g i+ 1 2 ,j+ 1 2 -g i+ 1 2 ,j-1 2 )f j -f i+1 f N g i+ 1 2 ,N - 1 2
).

If we bring together the terms f i+1 we get

K i+ 1 2 = ∆ε i+ 1 2 ∆ε i+ 1 2 (( 1 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2 + N -1 j=1 (g i+ 1 2 ,j+ 1 2 -g i+ 1 2 ,j-1 2 )f j -f N g i+ 1 2 ,N -1 2 )f i+1 - N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2 f i ∆ε i+ 1 2 ).
We denote Remark 20. One could replace the exact equilibrium M i in the definition of the scheme by an approximate equilibrium state M i with same mass and energy and the resulting scheme is still conservative, positive but the equilibrium state is now M i . This could be useful since this scheme requires the knowledge of the equilibrium state and normally this can be done only by solving a non-linear equation with an iterative method.

A i+ 1 2 = ∆ε i+ 1 2 ∆ε i+ 1 2   1 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2 + N -1 j=1 (g i+ 1 2 ,j+ 1 2 -g i+ 1 2 ,j-1 2 )f j -f N g i+ 1 2 ,N -1 2   , (3.

60) and

B i+ 1 2 = 1 ∆ε i+ 1 2 N -1 j=1 g i+ 1 2 ,j+ 1 2 f j+1 ∆ε j+ 1 2 , ( 3
Remark 21. If we stand

f i+ 1 2 = f i instead of f i+ 1 2
= f i+1 in the numerical flux and if we take

∆ε i+ 1 2 = (∆M ) i+ 1 2 βM i
, the Maxwellians are also preserved but we can show that S 3 is positive if and only if the energy grid is uniform.

Remark 22. On a uniform grid all the terms ∆ε i+ 1 2 are equal and up to a multiplicative constant the flux reduces to

K i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 (f i+1 f j -f i f j+1 ).
The scheme is then nothing else than the scheme provided by Berezin and Pekker [START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF] which is also an entropic scheme as shown in [START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF]. can be achieved in only O(N ) operations, as explained in [START_REF] Buet | Conservative and entropy schemes for the isotropic Fokker-Planck-Landau equation[END_REF][START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF].

Remark 24. The second idea of Larsen et al. in [START_REF] Larsen | Discretization methods for one-dimensional Fokker-Planck operators[END_REF] to preserve Maxwellians and applied here to the isotropic FPL equation is to write

∂f ∂ε = M ∂ ∂ε ( f M ) -βf . The isotropic FPL equation reads ∞ 0 ∂f (ε) ∂t φ(ε) √ εdε = - 1 2 ∞ 0 ∞ 0 ( ∂φ(ε) ∂ε - ∂φ(ε ) ∂ε )g(ε, ε ) f (ε )M (ε) ∂ ∂ε ( f (ε) M (ε) ) -f (ε)M (ε ) ∂ ∂ε ( f (ε ) M (ε ) ) dε dε.
Proceeding as above leads to a scheme which is indeed conservative in mass and energy and preserves the Maxwellians. But at this time, we cannot show the positivity. Nevertheless it could be possible to derive a positive scheme since there are sufficiently degrees of freedom in the discretization, namely the factor f i+ 1 2 and M i+ 1 2 .

Entropy decaying scheme (scheme S 4 )

The present approach consists in deriving a discretization from the "Log" form (2.4) because this helps to check easily the main properties of the operator, conservation and H-theorem which are given without proof. This approach was initiated by Degond and Lucquin in [START_REF] Degond | An entropy scheme for the Fokker-Planck collision operator of plasma kinetic theory[END_REF] for the 3-D Landau equation (for the implementation see [START_REF] Buet | Fast algorithms for numerical, conservative, and entropy approximations of the Fokker-Planck-Landau equation[END_REF][START_REF] Lemou | Multipole expansions for the Fokker-Planck-Landau operator[END_REF]), used for the 2-D axisymetric Landau equation [START_REF] Lemou | Multipole expansions for the Fokker-Planck-Landau operator[END_REF], and applied in the isotropic case in [START_REF] Buet | Numerical analysis of conservatives and entropy schemes for the Fokker-Planck-Landau equation[END_REF][START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF]. We present a new version of the algorithm based on the "Log" form.

From the previous section we retain the approximation of

ε 0 0 ∂f ∂t φ √ εdε as N i=1 c i df i dt φ i .
We now turn to the discretization of the right-hand side of (2.4) (r.h.s.) = -1 2

N -1 i=1 N -1 j=1 ε i+1 ε i ε j+1 ε j ( ∂φ(ε) ∂ε - ∂φ(ε ) ∂ε )g(ε, ε )f (ε )f (ε)( ∂ log f (ε) ∂ε - ∂ log f (ε ) ∂ε )dεdε .
(3.62) For each integrals of (3.62) we use again a midpoint quadrature formula, the ε-derivative are approximated by centered finite difference operator, thus (r.h.s.) = -1 2 where the k i+ 1 2 are taken as an arithmetic mean of f i and f i+1 . This yields a discrete model for which it cannot be proved the f remains positive as it must be. More recently, Cordier and Buet [START_REF] Buet | Conservative and entropy schemes for the isotropic Fokker-Planck-Landau equation[END_REF] have chosen a second-order approximation as the harmonic average; that is (2f i f i+1 )/(f i + f i+1 ). They proved the existence of a unique, positive and entropy solution for the semi-discretized FPL equation. Unfortunately, it turns out if initial data is not strictly positive. Furthermore, if the distribution function is zero at a point all numerical fluxes vanish. In the case of a uniform grid, in [START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF], and from an idea of [START_REF] Buet | Numerical solution of an ionic Fokker-Planck equation with electronic temperature[END_REF], the authors have considered the formula

N -1 i=1 N -1 j=1 ∆φ i+ 1 2 ∆ε i+ 1 2 - ∆φ j+ 1 2 ∆ε j+ 1 2 g i+ 1 2 ,j+ 1 2 k i+ 1 2 ,j+ 1 2 (∆ log f ) i+ 1 2 ∆ε i+ 1 2 - (∆ log f ) j+ 1 2 ∆ε j+ 1 2 ∆ε i+ 1 2 ∆ε j+ 1 2 , ( 3
k i+ 1 2 ,j+ 1 2 = f i (∆f ) j+ 1 2 -f j (∆f ) i+ 1 2 (∆ log f ) j+ 1 2 -(∆ log f ) i+ 1 2 if (∆ log f ) i+ 1 2 = (∆ log f ) j+ 1 2 ,
and

k i+ 1 2 ,j+ 1 2 = f i f j when (∆ log f ) i+ 1 2 = (∆ log f ) j+ 1 2 .
For a uniform grid and only in this case the above expression can be simplified into

k i+ 1 2 ,j+ 1 2 = f i f j+1 -f j f i+1 log(f j+1 f i ) -log(f i+1 f j ) ,
with such assumption they recover the scheme proposed in [START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF]. Following this approach we take As in the previous sections the system of ordinary equation reads

df i dt = Q S 4 i , 1 ≤ i ≤ N (3.64)
where

Q S 4 1 = K 3 2 /c 1 , Q S 4 i = (K i+ 1 2 -K i- 1 2
)/c i for 2 ≤ i ≤ N -1 and

Q S 4 N = -K N -1 2
/c N with the numerical flux

K i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 ( f i+1 + f i 2 ) 1-a i+ 1 2 ( f j+1 + f j 2 )
1-a j+ 1 Remark 25. On a uniform grid the flux reduces to

K i+ 1 2 = N -1 j=1 g i+ 1 2 ,j+ 1 2 (f i+1 f j -f i f j+1 ),
thus the scheme reduces to the scheme of Berezin and Pekker [START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF] which is also an entropic scheme as shown in [START_REF] Buet | Numerical analysis of the isotropic Fokker-Planck-Landau equation[END_REF].

In this work, we have analyzed the discretization of the isotropic FPL equation in order to obtain positive, conservative and Maxwellian state preserving schemes. Despite the lack of existence and uniqueness results for the schemes presented here, we have clarified the derivation as well as the properties of the Chang and Cooper method in this nonlinear case. We have clearly explained how bad boundary conditions could lead to negative or non conservative schemes. As we have seen, the Chang and Cooper method is a quite complicated scheme for a non linear equation such as the FPL equation. We have also shown that the Epperlein scheme could not guarantee energy conservation and Maxwellian preserving, and the scheme S 2 (developed by Langdon in the case of Coulombian potentials [START_REF] Decoster | [END_REF][START_REF] Langdon | Conservative differencing of the electron Fokker-Planck transport equation[END_REF]), is not a positive scheme. Only S 1 , the new scheme we developed, fulfills all the requirements. Let us recall that on infinite grids these three schemes are identical. The differences stem from the truncation of the domain of computation and the choice of the boundary conditions. In conclusion a direct discretization of the special form (2.15) and (2.16) of the drift and diffusion coefficients is not a guarantee to ensure energy conservation as claimed by Kingham and Bell in [START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF].

We have also presented and analyzed in detail two simpler schemes, namely S 3 and S 4 which are not based on the Chang and Cooper method. However, they are as efficient as S 1 . The scheme S 4 based on the "Log" form of the operator is also an entropic scheme. Let us also mention that the Logarithm of the distribution function would not be a problem for an implicit time discretization since the solution of the FPL equation would be positive for all time t > 0 as soon as the initial data is non-negative. We can summarize the properties of the above schemes as follow:

• Chang and Cooper S 1 : positive, conserves the energy and the Maxwellians.

• Chang and Cooper S 2 : non-positive, conserves the energy and the Maxwellians.

• Chang and Cooper by Epperlein: positive, does not conserve the energy and the Maxwellians.

• S 3 : positive, conserves the energy and the Maxwellians.

R

  A G N É T I S M E D I R E C T I O N D E S S Y S T È M E S D ' I N F O R M A T I O N C E A / S A C L AY 9 1 1 9 1 G I F -S U R -Y V E T T E C E D E X F R A N C E

1 2(positive or negative) if δ i+ 1 2 exists we get 0 ≤ δ i+ 1 2 ≤ 1 . 2 '

 12212 Therefore, since the f i 's are positive the f i+ 1 s remain positive. Now we clarify the expression of δ i+ 1 2

1 2

 1 is given by the Chang and Cooper average (3.11)-(3.12).

3. 1 . 2

 12 The scheme S 2 (Langdon's scheme in the case of a Coulombian potential) Another option consists in integrating the right-hand side of (3.6) up to ε N +1 , assuming that f N +1 = 0 and constraining the flux at the last point ε N + 1 2

1 2

 1 is given by the Chang and Cooper average (3.23)-(3.24).

1 2is undefined. As the f i+ 1 2 '

 12

  .53) By summing the (N -3) relations (3.52) together with (3.53) we easily get the equation at i = 1. This ends the proof.

2 's and B i+ 1 2 ' 1 2≥ 0 . 2 and B i+ 1 2 2 and B i+ 1 2

 22102122 s are positive and bounded. As a result (see Remark 6) we have D i+ Moreover u and v are positive functions therefore A i+ 1 are positive too. Now we have to show that A i+ 1

1 2and B i+ 1 2 2 ) 2 ) 2 and B i+ 1 2

 112221 remain bounded. Now we assume that | α i+ 1 2 |≤ 1 then u(α i+ 1 and v(α i+ 1 are bounded and A i+ 1 too since D i+ 1 2

  Now we show that D i+ 1 2 behaves like 1/Y . First, we saw that the drift coefficients are always positive and the diffusion coefficients (namely D i+ 1 )) have all the same sign. Consequently the Chang and Cooper coefficients always verify: 0

3. 1 . 8

 18 Discussion: some remarks and about Epperlein's version of the Chang and Cooper method.

.61) thus we have K i+ 1 2 = A i+ 1 2 f i+1 -B i+ 1 22 ,j+ 1 2 2 ≤ 1 (g i+ 1 2 ,j+ 1 2 2 ≤ ∆ε N - 1 2 2 . Now, ∆ε i+ 1 2 /∆ε i+ 1 2 ∈ 2 , B i+ 1 2

 221121112122221 f i .We can summarize the properties of this scheme in the following proposition Proposition 9. If g i+ 1 is an increasing sequence and ∆ε i+ 1 ∆ε N -1 2 the scheme S 3 is conservative in mass and energy. Moreover it is positive and preserves the equilibrium state when is reached. Proof. Mass and energy conservation are easily checked by setting φ = 1 and φ = ε in the weak symmetrized form of the discrete FPL equation (3.58). Using the definition (3.57) of ∆ε i+ is an increasing sequence and that ∆ε i+ 1 leads to the positivity of A i+ 1 [0, 1] therefore A i+ 1 are bounded since, due to mass conservation, we have f i ≤ ρ min j c j . According to Lemma 1, f i cannot vanish in finite time. That completes the proof.

.63) with g i+ 1 2 ,j+ 1 2 = g(ε i+ 1 2 , ε j+ 1 2 ) and the terms k i+ 1 2 ,j+ 1 2 2 f j+ 1 22 ,j+ 1 2 are of the form k i+ 1 2 ,j+ 1 2 = k i+ 1 2 k j+ 1 2

 1222112121212 stand for a second-order approximation of f i+ 1 at the center of the interval [ε i , ε i+1 ] × [ε j , ε j+1 ]. In the paper of Berezin et al.[START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF], the terms k i+ 1

k i+ 1 2 ,j+ 1 2 = l i+ 1 2 2 and for consistency we take l i+ 1 2 = ( f i + f i+1 2 ) 1 -a i+ 1 2 . 1 2 ,j+ 1 2

 1212221211 By the relation (x -y)(log x -log y) ≥ 0 we assure that k i+ remains positive. Thus we obtain the following weak formulation for the semi-discretized model

2 =i+ 1 2 -1 and Θ l i+ 1 2 =i+ 1 2 = N - 1 j=1 g i+ 1 2 ,j+ 1 2 (Θ r i+ 1 2 f i+1 Θ l j+ 1 2 f j -Θ l i+ 1 2 f i Θ r j+ 1 2 2 = A i+ 1 2 f i+1 -B i+ 1 2 f i . Proposition 10 .i+ 1 2 = A i+ 1 2 f i+1 -B i+ 1 2 f i , where A i+ 1 2 and B i+ 1 2 2 = 1 2- 1 ≤ 2 = 1 2- 1 ≤ 1 thus A i+ 1 2 and B i+ 1 2

 222211222212211022121211211121 (2f i+1 /(f i + f i+1 )) a (2f i /(f i + f i+1 ))The scheme S 4 is conservative in mass and energy and preserves the positivity and the Maxwellian equilibrium when it is reached. Moreover it is an entropic scheme.Proof. The conservation of mass and energy are directly checked on (3.63) by taking φ = 1 or φ = ε. Let M the Maxwellian equilibrium solution. The proof is straightforward by substituting the distribution function M also in(3.63). If we define the entropy byH(f ) = N i=1 c i f i log f i , thus using (3.63) with φ = log f , one has dH(f ) dt = N i=1 c i df i dt log f i ≤ 0,and the equality occurs when f = exp(-βε), that is when f is a Maxwellian. In other terms we have a H-theorem for this scheme. The numerical flux associated to scheme S 4 is K are respectively given by (3.65) and (3.66). If{f i } 1≤i≤N is positive it's straightforward to show that 0 ≤ Θ l i+ 1 (2f i /(f i+1 + f i )) a i+ 1 and 0 ≤ Θ r i+ 1 (2f i+1 /(f i+1 + f i )) a i+are positive and upper bounded since, due to mass conservation, we have f i ≤ ρ min j c j . Thus according to Lemma 1 the solution cannot vanish in finite time and consequently it is positive. That ends the proof.

  .32) Note that d(f ) doesn't depend on (i +1 2 ). Thus, the Chang and Cooper relation is reduced to a simple scalar non-linear equation. To find δ i+ 1 2 all we need to do is to solve equation (3.32) and that can be achieved by a Newton method in only O(N ) operations at each step. Moreover if the grid is uniform all the α i+ 1 2 and δ i+ 1 2 are equal and one can check directly on (3.6) that the scheme reduces to the scheme of Berezin and Pekker[START_REF] Yu | Conservative finite difference schemes for the Fokker-Planck equation not violating the law of an increasing entropy[END_REF]

  .38) By construction(3.36) and(3.38) are consistent approximations of drift E(f ) and diffusion D(f ) when their integration domain is reduced to a bounded domain [0, E] in the variable ε (see equations (2.13) and (2.14)). Now, the collision terms computed by the scheme S 2 read in the case of Coulombian interactions

  .44)

	Remark 11. In the Coulombian case, thanks to (3.35), (3.39), (3.37), (3.41) and (3.21) we
	can write

  j |< 1 and H is contracting. Now if we admit that β is sufficiently small and ∆ε/∆ε i+ 1

	2	close to the unity such that | g(α i+ 1 2	2 )∆ε/∆ε i+ 1	|< 1 then H is still contracting.

f (ε) ,(2.13)

Acknowledgments

We wish to thank A. Decoster for many useful discussions during the course of this work.

• S 4 : positive, conserves the energy and the Maxwellians and is entropic.

We hope that some of the remarks would help people to design an implicit time discretization both positive and conservative. Concerning the implementation of implicit time discretization, our conviction is that a fixed point method based on the form of the fluxes

at the previous iteration is sufficient. At each iteration there is no energy conservation (see also [START_REF] Epperlein | Implicit and conservative difference schemes for the Fokker-Planck equation[END_REF]) but the iteration is positive, and by making the fixed point converge to zero machine size to enforce energy conservation would be efficient enough. This is the method chosen in [START_REF] Buet | A Conservative and Entropy Scheme for a Simplified Model of Granular Media[END_REF] in the context of granular media for the same type of equation and by Kingham and Bell in the code IMPACT for the isotropic FPL equation [START_REF] Kingham | An implicit Vlasov-Fokker-Planck code to model nonlocal electron transport in 2-D with magnetic fields[END_REF]. But people could also follow the alternative ideas developped by Lemou and Mieussens in [START_REF] Lemou | Implicit schemes for the Fokker-Planck-Landau equation[END_REF] for implicit time discretization of the 3-D or isotropic Landau equation.

LA DIRECTION DES SYSTEMES D'INFORMATION

CEA / SACLAY 91191 GIF-SUR-YVETTE CEDEX FRANCE