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-Landau 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.


Ngoc Q K Duong

Alexey Ozerov

Louis Chevallier

Joël Sirot

AN INTERACTIVE AUDIO SOURCE SEPARATION FRAMEWORK BASED ON NON-NEGATIVE MATRIX FACTORIZATION

Keywords: Interactive audio source separation, nonnegative matrix factorization, uncertainty-based learning, timefrequency annotation, user feedback

Though audio source separation offers a wide range of applications in audio enhancement and post-production, its performance has yet to reach the satisfactory especially for single-channel mixtures with limited training data. In this paper we present a novel interactive source separation framework that allows end-users to provide feedback at each separation step so as to gradually improve the result. For this purpose, a prototype graphical user interface (GUI) is developed to help users annotating time-frequency regions where a source can be labeled as either active, inactive, or well-separated within the displayed spectrogram. This user feedback information, which is partially new with respect to the state-of-the-art annotations, is then taken into account in a proposed uncertainty-based learning algorithm to constraint the source estimates in next separation step. The considered framework is based on non-negative matrix factorization and is shown to be effective even without using any isolated training data.

INTRODUCTION

Audio source separation still remains a very hot research topic even there has been a surge of research over the past years [START_REF] Vincent | The Signal Separation Campaign (2007-2010): Achievements and remaining challenges[END_REF]. In fact, separation performance that can be achieved by state-ofthe-art algorithms is far from satisfactory in certain scenarios such as under-determined mixtures of reverberated sources. In singlechannel case, where spatial information about the sources cannot be exploited, the problem becomes even harder and usually relevant training data is needed to first learn the spectral characteristics of individual sources. Such a class of supervised algorithms is mostly based on non-negative matrix factorization (NMF) [START_REF] Virtanen | Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria[END_REF][START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF], or its probabilistic formulation known as probabilistic latent component analysis (PLCA) [START_REF] Raj | Latent variable decomposition of spectrograms for single channel speaker separation[END_REF][START_REF] Smaragdis | Supervised and semi-supervised separation of sounds from single-channel mixtures[END_REF]. However, when training examples are unavailable (or not representative enough) these methods can not be applied without other prior information about the sources. Examples of such prior information include "hummed" sounds that mimic the ones in the mixture [START_REF] Smaragdis | Separation by humming: User-guided sound extraction from monophonic mixtures[END_REF], or text transcriptions of the corresponding speech sources [START_REF] Magoarou | Textinformed audio source separation using nonnegative matrix partial co-factorization[END_REF].

So called user-guided approaches based on NMF have been proposed recently. These approaches allow end-user to manually annotate information about activity of each sound source, e.g. if it is active or not, in time [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF] or time-frequency (T-F) [START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF] domains. The annotated source activity information is then used, instead of training data, to guide the separation process. Another user-guided approach based on independent vector analysis, which allows user to tune temporal power variations of sources, was also presented in [START_REF] Ono | User-guided independent vector analysis with source activity tuning[END_REF]. Though prior work has shown the effectiveness of such algorithms, the results are still far from perfect especially in mixing scenarios with strong overlapping sources. A more advanced interactive framework based on PLCA was proposed by Bryan et al. [START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF] whereby overall separation process comprises several interactive separation steps (if needed). At each step, end-user can perform T-F annotation on the spectrogram displays of intermediate separation results, in addition to the annotation on the spectrogram of the mixture itself. Note that in contrast to [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF][START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF], where the annotations are specified only once, interactive approach in [START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF] allows user feedback at each separation step so as to gradually improve the result by correcting remaining errors.

Motivated by the effectiveness of both the user-guided algorithms [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF][START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF] and the interactive strategy [START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF] 1 , we present in this paper an interactive source separation approach based on a novel NMF formulation. The proposed approach can efficiently handle user feedback information at each separation step, and is robust to errors in the annotations thank to the derived uncertainty-based learning algorithm for parameter estimates. In addition to the T-F annotation about the source activity considered in the existing algorithms, we introduce a new type of annotation about the quality of the separated sources: user can validate if a source is well-separated in certain T-F regions. These well-separated regions can then be effectively exploited, as sort of training data, in our formulation to better constrain parameter estimates in next separation step.

The rest of the paper is organized as follows. In Section 2 we summarize the NMF model for source separation and a baseline parameter estimate exploiting temporal annotation as in [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF]. This baseline can serve as the first separation step in our interactive framework. Main contributions of the paper are described in Section 3 where the T-F annotations via a GUI are introduced followed by the proposed optimization algorithm and a global description of the interactive system. We conduct experiments to validate the effectiveness of the proposed approach in Section 4. Finally we conclude in Section 5.

NMF MODEL AND BASELINE PARAMETER UPDATE

Model

Let x f n and s j,f n be the short-time Fourier transform (STFT) coefficients of the observed single-channel mixture signal and the con-tribution of j-th source signal, respectively, at frequency bin f and time frame n. The mixing model writes

x f n = J j=1 s j,f n , ( 1
)
where J is the total number of sources, f = 1, . .

Baseline parameter update exploiting temporal annotation

NMF parameters can be estimated in the maximum likelihood (ML) sense, which is equivalent to minimizing the following cost function

[3] C1(θ) = f n dIS (v f n |v f n ) , (2)
where θ = {W, H} being the NMF parameters, and dIS (•|•) denotes Itakura-Saito (IS) divergence. The well-known multiplicative update (MU) rules for parameter estimation write:

H ← H W T (WH) .-2 V W T (WH) .-1 (3) W ← W (WH) .-2 V H T (WH) .-1 H T (4)
where denotes the Hadamard entrywise product, A .p being the matrix with entries [A] p ij , and the division is entrywise. Considering the baseline user-guided approach [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF], parameter estimation is guided by temporal annotation where user is asked to indicate time segments along the mixture where each source is active. This temporal annotation is then reflected in the initialization of the activation matrix H: when j-th source is not annotated by "active" at time frame n, the n-th column of H (j) is initialized by 0 and it will not be changed via the multiplicative update. Note that this temporal annotation is much easier and faster than the detail time-frequency annotation considered in [START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF][START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF] as well as later part of the paper. Thus, with the absence of isolated training sounds, this baseline can serve as the first separation step in our interactive framework.

Source reconstruction

Given the estimated parameters θ = {W, H}, the source STFT coefficients are computed by Wiener filtering as ŝj,fn = vj,fn vfn

x f n , (5)
then the time domain source estimates are obtained via the inverse STFT.

PROPOSED INTERACTIVE SEPARATION ALGORITHM

We first present all types of the considered time-frequency annotations, which would bring benefit to the separation process, in Section 3.1. We then propose an uncertainty-based learning algorithm that allows to efficiently incorporate all annotated information to constrain parameter estimates in interactive separation steps in Section 3.2.

Time-frequency annotations

In order to support the annotation task, we developed a prototype user interface enabling user to select T-F regions (either in rectangle or polygon shapes) in the spectrogram display (either in linear or logarithmic frequency scale) of the mixture as well as the separated sources. This GUI is visualized in Fig. 1. Given a selected T-F region, and possibly by hearing the corresponding sound, end-user can assign to each source one of the following labels:

(i) source is active (default label),

(ii) source is inactive, (iii) source is well-separated.

Fig. 1. GUI for user annotation. Piano sound is labeled as wellseparated in the current selected polygon region (in blue) where phone-ring sound labeled as is inactive. User also annotated that speech is active alone (i.e. all the other sources are inactive) in the rectangle region (in yellow).

While annotation (i) and (ii) have been addressed in [START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF][START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF], annotation (iii) is new. We find that this type of information is interesting because first, the well-separated T-F regions should be constrained so as they will not be accidently damaged due to the changing of parameters in next separation steps. Second, and even more important, these regions can be used, as sort of input training data, to help better estimating the model parameters of the corresponding source.

Interactive parameter update exploiting T-F annotations

We derive an optimization algorithm addressing two challenges: (1) how to use annotated T-F regions to learn NMF models, and (2) how to handle errors in the annotation. While the former can be solved in the same way as in [START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF], the latter is more difficult. Also, the latter is crucial so as to account for the well-separated regions where artifacts and interferences still remain from the previous separation step2 . To fulfil this challenge, we exploit uncertainty-based learning [START_REF] Arberet | A tractable framework for estimating and combining spectral source models for audio source separation[END_REF][START_REF] Ozerov | Uncertainty-based learning of acoustic models from noisy data[END_REF] where source models can be learned from the source estimates, while taking into account separation errors described by their variances. The principle is summarized as follow. Within the above-mentioned Gaussian assumption, the posterior of j-th source writes [START_REF] Arberet | A tractable framework for estimating and combining spectral source models for audio source separation[END_REF] p(s j,f n |x f n ; θ) = Nc s j,f n ; ŝj,fn , vj,fn vj ,f n vfn ,

where vj ,f n = vfn -vj,fn and ŝj,fn is the reconstructed source from Wiener filtering [START_REF] Smaragdis | Supervised and semi-supervised separation of sounds from single-channel mixtures[END_REF], which can be then written as

ŝj,fn = s j,f n + b j,f n , (7)
where b j,f n ∼ Nc(0, vj,fn vj ,f n /v f n is source estimation error.

Estimating the NMF model for j-th source from ( 7) in the ML sense can be shown to be equivalent to minimizing the following cost [START_REF] Arberet | A tractable framework for estimating and combining spectral source models for audio source separation[END_REF] C2(θ) = j,f n dIS (ṽ j,f n |v j,f n + ẽj,fn ) ,

where ṽj,fn = (v j,f n /v f n ) 2 v f n and ẽj,fn = vj,fn vj ,f n /v f n are computed from model parameters estimated at the previous separation step and fixed at current estimation step. It can be seen that this parameter estimation takes into account both posterior mean ŝj,fn and posterior variance ẽj,fn characterizing source reconstruction error.

The cost ( 8) is now used to account for the annotated T-F regions in our framework. Combining (2) and ( 8), the overall proposed cost function to be optimized for interactive parameter estimates is

C3(θ) = f n dIS (v f n |v f n )+ j,f n λ j,f n dIS (ṽ j,f n |v j,f n + ẽj,fn ) , (9)
where λ j,f n , ṽj,fn and ẽj,fn are defined, for each source j at each time-frequency point (f, n), according to the T-F annotations as:

(λ j,f n , ṽj,fn , ẽj,fn ) = =          (a, 0, 0) if j is inactive, (b, v f n , 0) if j is active alone, c, vj,fn vfn 2 v f n , vj,fn vj ,f n vfn if j is well-separated,
(0, any, any) otherwise. [START_REF] Ono | User-guided independent vector analysis with source activity tuning[END_REF] where positive constants a, b, and c determine contributions of different annotation types, any means any value, and a source is active alone when all other sources are labeled as inactive. Note that all parameters vj,fn , vfn , and vj ,f n in [START_REF] Ono | User-guided independent vector analysis with source activity tuning[END_REF] are computed from the previous separation step and fixed in current interactive step. Further explanation for the quantities in ( 10) is as follow. When a source is labeled as inactive or active alone, its corresponding separation error variances ẽj,fn are set to zero, while ṽj,fn should be equal to either zero or the mixture power spectrum v f n , respectively. Besides, if a source is already well-separated in certain regions, its corresponding STFT coefficients should be constrained to keep values from the previous separation step and not to vary greatly. Let Λ (j) , V (j) and E (j) the F × N matrices characterizing the T-F annotations for j-th source with, respectively, λ j,f n , ṽj,fn and ẽj,fn as entries. Using general principle from [START_REF] Févotte | Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis[END_REF] we derived MU rules to minimize [START_REF] Lefèvre | Semi-supervised NMF with time-frequency annotations for singlechannel source separation[END_REF] and the resulting parameter updates are as follows:

H (j) ← H (j) W T (j) WH .-2 V + Λ (j) W (j) H (j) + E (j) .-2 V (j) W T (j) WH .-1 + Λ (j) W (j) H (j) + E (j) .-1 (11)
W (j) ← W (j) WH .-2 V + Λ (j) W (j) H (j) + E (j) .-2 V (j) H T (j) WH .-1 + Λ (j) W (j) H (j) + E (j) .-1 H T (j) (12)

Overall system description

Following the baseline separation step summarized in Section 2.2 and the proposed interactive refinement steps with user's T-F annotations introduced in Section 3.2, this section summarizes the global workflow of the overall system. For clarity, all steps in the proposed interactive framework are presented in Algorithm 1.

Note that all T-F annotations in interactive separation steps are kept in the memory and will be reused, together with new annotations, in next separation step.

EXPERIMENTS

We first evaluated the source separation performance of the proposed approach over four single-channel mixtures of different male and female speeches (extracted from TIMIT database) with different background sounds (piano chords, drums, and phone ring). The mixtures were 23 second long, sampled at 16 kHz, and articially mixed at 0 dB SNR. In each mixture speeches (either male or female voices) and background sounds appear alone during about three seconds so that temporal annotation was well-exploited to guide the parameter estimates in the first separation step.

We compared separation result of the derived interactive approach after the first separation step (named Int-SS-1) with that obtained after the second separation step when the proposed annotation about the quality of separated sources, i.e. well-separated T-F regions, is either not used (named Int-SS-2) or used (named Int-SS-3). The time-frequency annotation was performed by one of the author who has experience on sound processing. Note that Int-SS-1 is equivalent to the user-guided algorithm [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF] where only temporal annotation was used, and Int-SS-2 is comparable to the state- of-the-art interactive method [START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF] though the former is formulated based on NMF while the later uses PLCA. For parameter estimation, the number of MU iterations was set to 100 for all separation steps, Kj = 20, j = 1, 2, and the NMF parameters W (j) , H (j) were re-initialized, in the second separation step, in the same way as they were initialized in the first step. We observed that this reinitialization brings slightly better result than initializing directly by the parameters obtained after MU iterations in the first separation step. The trade-of parameters a, b, and c are set manually between 1 and 10.

Approach SDR SIR SAR Int-SS-1 (comparable to [START_REF] Ozerov | Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation[END_REF]) 8.9 16.0 13.1 Int-SS-2 (comparable to [START_REF] Bryan | Interactive refinement of supervised and semi-supervised sound source separation estimates[END_REF] Separation performance was evaluated using the signal-todistortion ratio (SDR) criterion measuring overall distortion, as well as the signal-to-interference ratio (SIR), and signal-to-artifact ratio (SAR) criteria [START_REF] Vincent | The 2008 Signal Separation Evaluation Campaign: A community-based approach to large-scale evaluation[END_REF], measured in dB, averaged over all sources, and shown in Table 1. As expected, performing source separation interactively even with only one user-feedback step (Int-SS-2 and Int-SS-3) significantly improves the result, i.e. in terms of SDR and SIR, compared to the baseline (Int-SS-1). This is because some interferences remained in separated sources from the first separation step, which were annotated by user, were successfully removed in the second step. Besides, it is worth noticing that Int-SS-3 improves the SDR, SIR, and SAR by 0.4 dB, 0.6 dB, and 0.2 dB, respectively, compared to Int-SS-2. This result confirms the effectiveness of our proposed well-separated annotation. In addition to the above synthetic mixtures, we have also tested the proposed approach on a twelve minute real-world sound track of Beverly Hills drama series for the separation of dialogs (speeches from main actors and actresses) and background sounds (music, environmental noise, footstep, etc). Due to the lack of groundtruth signals for objective evaluation, informer listening has confirmed the quality of separated sources after two interactive steps with about ninety minutes of manual annotation.

CONCLUSION

In this paper, we have presented a novel interactive source separation framework based on NMF formulation. The proposed approach efficiently exploits a new type of user annotation about the quality of the estimated sources, in addition to the existing T-F annotations about the contributions of the sources to the mixture spectrogram, to further constrain the parameter estimates in interactive separation steps. Preliminary experiments with both simulated and real-world mixtures confirm the effectiveness of the derived algorithm. Future research would be devoted to learn an automatic or semi-automatic annotation algorithm so as to reduce manual effort and fasten annotation time.

  . , F and n = 1, . . . , N . Defining the power spectrogram of the mixture by v f n = |x f n | 2 , NMF aims at approximately factorizing the F × N matrix V = {v f n } f n as V ≈ V = WH, where W and H are nonnegative matrices of size (F × K) and (K × N ), respectively. Assuming that the source STFT coefficients follow zero-mean Gaussian distribution s j,f n ∼ Nc(0, vj,fn ) where vj,fn = [W (j) H (j) ] f n , W (j) and H (j) are matrices of size (F × Kj) and (Kj × N ), respectively, modeling the contribution of the j-th source. With linear model vfn = j vj,fn [3], we can write W = [W (1) , . . . , W (J) ] and H = [H T (1) , . . . , H T (J) ] T .

  Algorithm 1 Global workflow of the proposed interactive source separation algorithm.1. User listens to the mixture and roughly annotate temporalsegments where each source is active 2. Randomly initialize NMF parameters {W (j) , H (j) }j, then set n-th column of H (j) to 0 if j-th source is inactive at time frame n (see Section 2.2)

	3. Perform first source separation step: parameters are alternat-
	ingly updated by (3) and (4) until convergence, then sources
	are reconstructed by Wiener filtering (see Section 2.3)
	4. User evaluates the separation performance. Stop separation
	process when the result is satisfied, otherwise:
	5. User performs the detailed T-F annotations (see Section 3.1)
	6. Construct the F × N nonnegative matrices Λ (j) , V (j) and
	E (j) for all source j from all the T-F annotations as (10)
	7. Initialize NMF parameters similarly to step 2
	8. Perform interactive source separation step: parameters are al-
	ternatingly updated by (11) and (12) until convergence, then
	sources are reconstructed by Wiener filtering
	9. Go to step 4.

Table 1 .

 1 Average source separation performance.

	)	9.8	16.8 13.1
	Int-SS-3 (proposed)	10.2 17.4 13.3

This approach obtained better objective separation result compared to other submissions in "Professionally produced music recordings" task of the Signal Separation Evaluation Campaign (SiSEC2013): http://www.onn.nii.ac.jp/sisec13/evaluation_ result/MUS/testMUS2013.htm

Though well-separated regions are validated by end-user via listening, they are usually not error-free. Indeed, artifacts and interferences still exist in these regions, but they are simply masked by the target source.


Abbas Bradai
email: bradai@labri.fr

Toufik Ahmed

Samir Medjiah
email: medjiah@labri.fr

QoE Assessment for SVC Live Distribution and Adaptation: ENVISION case study

Keywords: Quality of Experience, Streaming, P2P, Smoothing, SVC I

Scalable video coding has drawn great interest in content delivery in many multimedia services thanks to its capability to handle terminal heterogeneity and network conditions variation. In our previous work, and under the umbrella of ENVISION, we have proposed a playout smoothing mechanism to ensure the uniform delivery of the layered stream, by reducing the quality changes that the stream undergoes when adapting to changing network conditions. In this paper we study the resulting video quality, from the final user perception under different network conditions of loss and delays. For that we have adopted the Double Stimulus Impairment Scale (DSIS) method. The results show that the Mean Opinion Score for the smoothed video clips was higher under different network configuration. This confirms the effectiveness of the proposed smoothing mechanism.

INTRODUCTION

The popularity of digital video content together with the increasing number of connected end-users through heterogeneous networks motivates the utilization of Peer-to-Peer (P2P) networks. In these networks, the ability to provide a high level of user experience is a challenging task. This task is supported by effective evaluation of service quality. In order to realize such a quality evaluation service, an assessment model is required to bridge relevant influencing factors with end users' experience.

The typical approaches to video quality assessment are mostly based on Quality of Service (QoS). QoS is referred as a performance measure from the perspective of system performance or reliability. However, QoS metrics do not take into account human users' perception, therefore there is a lack of correlation between QoS and users' opinions. The concept of Quality of Experience (QoE) has been introduced [START_REF] Jain | Quality of experience[END_REF] to fill this gap between service performance and user experience. QoE is commonly defined as the overall acceptability of a service as subjective perceived by end users. Using models based on the concept of QoE, it is possible to measure the extent to which a user-centric service (e.g. video content distribution) achieves its objectives.

In our previous works [START_REF] Bradai | On the optimal scheduling in pull-based real-time p2p streaming systems: layered and non-layered streaming[END_REF][3] [START_REF] Bradai | An efficient playout smoothing mechanism for layered streaming in P2P networks[END_REF], under the umbrella of the European project ENVISION [5], we discussed several complexities involved in layered streaming over P2P networks due to bandwidth fluctuation and peer's unreliability. In this case, the correct decision regarding the selection and the scheduling of the layers is crucial: how many layers to request, in which order and from which peer? In [START_REF] Bradai | An efficient playout smoothing mechanism for layered streaming in P2P networks[END_REF], we proposed amplitude and frequency reduction (Figure 1) to ensure the smooth quality of the stream and bandwidth allocation schemes for layered streaming systems for P2P networks in [START_REF] Bradai | Differenciated bandwidth allocation in P2P layered streaming[END_REF]. The amplitude reduction mechanism aims to reduce the transition from the higher layer to very lower layer (amplitude reduction) and vice versa. The video quality assessment for amplitude reduction can be performed by examining the optimal trade-off between the temporal and SNR scalability. The video basement method such as Double Stimulus Impairment Scale (DSIS) [START_REF] Larabi | Quality Assessment Approaches[END_REF] can be used with different frame rate and quantisation parameter but having fixed spatial resolution. The experiment can be performed for different content type such as content having high motion, low motion etc.

Similarly, the subjective evaluation can be performed for frequency reduction mechanism. The goal of frequency reduction is to reduce the number of changes from one layer to another that affects user's perception of stream as described in [START_REF] Bradai | An efficient playout smoothing mechanism for layered streaming in P2P networks[END_REF].

II. STATE OF THE ART

Human perception of different temporal resolutions has been investigated for a relatively long time. Especially, much work has been done to find the minimum acceptable frame rates of video stimuli for various tasks such as target tracking, target detection/recognition, lip reading, orientation judgment, etc. [START_REF] Chen | Review of low frame rate effects on human performance[END_REF]. As for the scenario of video consumption, various factors affect the perceived quality of different frame rates, e.g. content type, viewing condition, display type, and observer characteristics.

For a certain target bit rate, once the spatial resolution is fixed, the temporal resolution can be increased only at the cost of decreased frame quality. The video assessment is performed in order to examine the optimal trade-off between these two dimensions in terms of perceived quality of the final video sequence. Traditionally, it is believed that a high frame rate is more important for content with fast motion than a high frame quality, which is supported by [START_REF] Yadavalli | Frame rate preferences in low bit rate video[END_REF]. In [START_REF] Yadavalli | Frame rate preferences in low bit rate video[END_REF], subjective experiments were conducted for video sequences encoded by using three different codec's (i.e. the Sorenson codec 2.1, H.263+ and a wavelet-based codec) for eight content types. For a fixed resolution of 352 × 240 pixels, three frame rates (10, 15 and 30 Hz) were considered. Overall, a frame rate of 15 Hz was most preferred across different coding conditions. However, content-dependence was observed, i.e. for content with slow (or fast) motion, preference of a frame rate of 10 Hz (or 30 Hz) was nearly as high as 15 Hz. Similarly, a double stimulus continuous quality scale (DSCQS) experiment in [WSV03] compared H.263+ sequences coded at three different frame rates (7.5, 15 and 30 Hz) with five quantisation parameter (QP) values, for a fixed spatial resolution of 320 × 192 pixels. The results showed that, for slow motion content, subjective quality degradation due to frame rate reduction was only minor. Different assessment methods have been recommended by the ITU [9], namely (a) The double-stimulus impairment scale (DSIS) [START_REF] Larabi | Quality Assessment Approaches[END_REF], (b) The double stimulus continuous quality-scale (DSCQS) [START_REF] Inácio | Quality user experience in advanced IP video services[END_REF], (c) Single-stimulus (SS) method [START_REF] Mantiuk | Comparison of four subjective methods for image quality assessment[END_REF], (d) Stimulus-comparison method [START_REF] Fiorucci | Objective and subjective quality assessment between JPEG XR with overlap and JPEG 2000[END_REF], (e) Single stimulus continuous quality evaluation (SSCQE) [START_REF] Liu | Continuous mobile video subjective quality assessment using gaming steering wheel[END_REF],

(f) Simultaneous double stimulus for continuous evaluation method (SDSCE) [START_REF] Bosc | Visual quality assessment of synthesized views in the context of 3D-TV. In 3D-TV System with Depth-Image-Based Rendering[END_REF].

In this paper we study the video quality from the final user perception for the SVC content smoothing in different network conditions of loss and delay, in ENVSION content distribution and adaptation framework. In this study we adopt the DSIS method, that will be described in detail in III.B. The rest of the paper is organized as follows: section III describes the evaluation environment, scenarios and methodology. Section VI presents the main obtained results and section V concludes this paper.

III. EVALUATION SCENARIOS AND METHODOLOGY

A. Evaluation environment and scenarios

The studied smoothing mechanisms was developed under the ENVISION framework. ENVISION proposes a crosslayer solution where the problem of supporting demanding services is solved cooperatively by service providers, ISPs, users and the applications themselves. To reach this goal, ENVISION relies on CINA [16] interface as a means for cooperation between applications and ISPs for realizing demanding multimedia services. -Use of CINA dynamic costs: The overlay content adaptation functions may respond to increased delay or decreased throughput caused in the event of congestion over an overlay link taking into account the CINA dynamic costs (Latency and LatencyLastHop policies).

-Increasing overlay versus network performance: The overlay topology construction algorithm may establish overlay connections using two different criteria: a) increase the overlay performance by connecting more peers to remote peers with good liveness (RoutingCost and Latency policies), or b) increase the network performance by reducing the number of connections to remote peers over paths with high CINA costs (RoutingCostLastHop and LatencyLastHop policies).

B. Evaluation method

Evaluation was made using the Double Stimulus Impairment Scale (DSIS) method, also referred as Degradation Category Rating (DCR). In this method videos are shown consequently in pairs. The first one is the reference, and expert is informed about it while the second one is impaired. After their playback, expert is asked to give his opinion on the second one, keeping in mind the first one. During different sessions (each session last up to half an hour) the assessors are shown series of pictures or sequences in random order and impairments covering all the essential combinations. It is important to note that the reference (unimpaired) picture is included in the pictures or sequences to be assessed. Finally, the mean score is calculated at the end of the session for each test condition and test picture.

In our evaluation, the original video is shown for 10 seconds, then the evaluated video for the same 10 seconds. This process was repeated for a total duration of 30 seconds of the video. The scores were expressed over 10, with 0 being the worst score and 10 being the best score.

IV. EVALUATION RESULTS

A. Performance enhancements at normal network conditions

At normal network conditions, i.e., no delays and no loss in the network, all the evaluated videos were rated above average for all the metrics used for the overlay construction except those based on Latency and LatencyLastHop. Figure 3 shows the obtained results. Also, we note for all the evaluated videos, smoothed videos achieved a higher MOS than videos without smoothing. In Figure 4, at the event of packet loss (10% loss) in the network, a good streaming experience (video quality) has been achieved for topologies constructed using "Random" and Latency metrics (> 6). However, we notice that random topology was more efficient in terms of MOS than topology built using Latency metric. Moreover, in this case, smoothed videos have not achieved a higher MOS than videos without smoothing.

C. Performance enhancements at the Event of Delays

At the event of delays in the network (delay of 500ms), the resulting streaming experience was rated above average (>5) in the case of different topologies. Figure 5 and Figure 6 show the obtained results. We note that videos from topologies built based on Latency and RoutingCost achieve the highest MOS (>6). Moreover, in general, smoothed videos are rated slightly higher than videos without smoothing. In this paper, we have presented a subjective quality measurement study of SVC Live content distribution and adaptation in ENVISION framework. The study was performed using the Double Stimulus Impairment Scale method by LaBRI students with some basic training in content adaptation. Overall, the Mean Opinion Score for the smoothed video clips was higher. In addition, the analysis of some of the low Mean Opinion Scores attributed to some clips has shown that even the slightest of quality fluctuations (a single missed frame) may have a big impact on the value of the Mean Opinion Scores (drop by two points). Thus, the effectiveness of the proposed smoothing mechanism.

Figure 1 :

 1 Figure 1: Variation in quality level (source [4])

Figure 2 :

 2 Figure 2: ENVISION overview In order to permit subjective user evaluation, the streaming application was run under different network and overlay scenarios in ENVISION testbed [17]: -Use of CINA static costs: The overlay topology construction algorithm may establish overlay connections randomly (Random policy), or taking into account the CINA static costs associated with the topology and remaining unchanged in the event of congestion (RoutingCost, RoutingCostLastHop policies).

Figure 3 :

 3 Figure 3: Mean Opinion Score -No delay, no loss B. Performance enhancements at the Event of Loss

Figure 4 :

 4 Figure 4: Mean Opinion Score -10% errors, no delay

Figure 5 :

 5 Figure 5: Mean Opinion Score -500ms delay, no errors

Figure 6 :

 6 Figure 6: Mean Opinion Score -1000ms delay, no errors

Acknowledgment

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) in the ENVISION project, grant agreement 248565. We want also to thank ENVISION partners: ULC, Orange, Alcatel-Lucent, LiveU and Telefonica to their contribution to the design of the FT testbed used in the evaluation scenario.