Christophe Buet
email: christophe.buet@cea.fr

Stéphane Cordier
email: stephane.cordier@univ-orleans.fr

Brigitte Lucquin-Desreux

Simona Mancini
email: smancini@ann.jussieu.fr

Éphane Cordier

Diffusion limits of the Lorentz model: Asymptotic preserving schemes

Keywords: Mathematics Subject Classification. 82C70, 35B40, 65N06 Hilbert expansion, diffusion limit

come

Introduction

In this article, we study the numerical schemes for a kinetic equation in the diffusive regime:

ε∂ t f +cosθ∂ x f = 1 ε L(f ). (1.1)
The problem is one-dimensional in the space variable x and bi-dimensional in the velocity variable v = (cos θ, sin θ). The unknown distribution function f = f (x, θ, t) is a function of position x (x ∈ R: in this article we shall not consider boundary conditions), of the velocity angle θ ∈ [-π, π]a n do ft i m et>0. The operator L(f ) is a linear collision operator of Lorentz type. Lorentz operators appear for example when considering elastic collisions of heavy particles (e.g. ions) against light ones (e.g. electrons); it is the first order term of the inter-species collision operator representing the collisions of the heavy particles on the light ones, when doing an asymptotic expansion in terms of the small mass ratio (see [START_REF] Degond | The asymptotics of collision operators for two species of particles of disparate masses[END_REF][START_REF] Lucquin-Desreux | Diffusion of electrons by multicharged ions[END_REF]). This operator does not depend on the energy variable, i.e. on the modulus of the velocity. It is defined in the Boltzmann case by (for simplicity, we drop the dependence with respect to x and t, since the collision operator only acts on θ):

L(f )(θ)= S 1 K(θ ′ -θ)[f (θ ′ ) -f (θ)] dθ ′ (1.2)
and in the Fokker-Planck case by:

L(f )(θ)=∂ 2 θθ f. (1.3)
We recall that the Boltzmann-Lorentz operator (1.2) converges (up to multiplication by the second order moment of the scattering cross-section K) to the Fokker-Planck-Lorentz one (1.3) when the cross-section concentrates, i.e. when the scattering angle during a collision, θ ′ -θ, is very small. This is the so-called "grazing collision limit" [START_REF] Buet | The grazing collision limit for the Boltzmann-Lorentz model[END_REF][START_REF] Cordier | Numerical approximation of the Vlasov-Fokker-Planck-Lorentz model[END_REF][START_REF] Lucquin-Desreux | A finite element approximation of grazing collisions[END_REF], which is also valid in the non-linear case [START_REF] Degond | The Fokker-Planck asymptotics of the Boltzmann collision operator in the Coulomb case[END_REF][START_REF] Desvillettes | On asymptotics of the Boltzmann equation when the collisions become grazing[END_REF]. It is well known in literature, that for ε ≪ 1 the solution of (1.1) converges to the solution of a diffusion problem, with respect to the space variable. The diffusion coefficient of this last equation may be computed by means of a Hilbert expansion method (see Sect. 2). Our goal is to derive a numerical scheme that is relevant for any value of ε. For simplicity, we will consider the Fokker-Planck-Lorentz collision operator defined by (1.3). Nevertheless, we remark that our results can be extended to the Boltzmann-Lorentz collision operator defined by (1.2).

This paper is inspired by a series of articles written by Jin and Levermore about diffusive limit of the isotropic Boltzmann-Lorentz operator, where the cross-section is constant with respect to the velocity variable. This fact does not allow to pass to the grazing collision limit, excluding thus the Fokker-Planck case. The authors consider discrete ordinate methods in the velocity variable, or equivalently in the cosine of the angle (the distribution function being isotropic in the others directions). More precisely, they determine the quadrature points in the integral with respect to velocity such that, in the diffusion limit, one recovers a heat equation with the correct diffusion coefficient. A first paper, see [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF], is devoted to the discretization in velocities, the distribution function depending continuously of the space and time variables. The authors construct quadrature sets corresponding to a small number of discrete velocities such that the diffusion coefficient and also the boundary conditions are compatible with the diffusive limit. In a second paper, see [START_REF] Jin | Fully-discrete numerical transfer in diffusive regimes[END_REF], they investigate the fully discrete case (both on the velocity and the space variables), the problem being stationary. This can be seen as the problem for one time step iteration using an implicit scheme. In another paper with Golse, see [START_REF] Golse | The convergence of numerical transfer schemes in diffusive regimes I: discrete-ordinate method[END_REF], they study the convergence of theses schemes.

Our aim is to extend this analysis to other elastic collision operators with an arbitrary cross-section (i.e. not necessarily constant) and for example in the case of the grazing collision limit, i.e. for the Fokker-Planck Lorentz operator (1.3).

There is a huge literature on related topics taking different names. We claim that the following keywords are used to qualify very close problematic although the goal and the methods used are different.

• Diffusive (a = 1) or hydrodynamic (a = 0) limit

ε a ∂ t f +cosθ∂ x f = 1 ε L(f ).
The aim is to derive numerical schemes which can be used either in rarefied (where ε is of order 1) or in dense regions (ε ≪ 1). Depending on the collision operator L, one expects the solution to converge towards an equilibrium state (typically a Maxwellian distribution) and that the conserved quantities are solution of the Euler equations (for a = 0) or Navier-Stokes equations (for a = 1), see [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF][START_REF] Klar | An asymptotic-induced scheme for non stationary transport equations in the diffusive limit[END_REF][START_REF] Prendergast | Numerical hydrodynamics for gas-kinetic theory[END_REF][START_REF] Prendergast | Numerical Navier-Stokes solutions from gas kinetic theory[END_REF].

A lot of works have been devoted to the coupling of rarefied and hydro-dynamical domains. For this purpose, numerical schemes compatible with the fluid limit, called Asymptotic preserving schemes, are an alternative to the matching of boundary conditions, these last being sometimes hard to design. In fact, these schemes can work uniformly with respect to the relaxation parameter (see [START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF] and its introduction for more details).

• Stiff source term for hyperbolic systems

A widely studied topic deals with hyperbolic systems of conservation laws (e.g. Euler equations) with stiff source terms (see [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF][START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]) due, for example, to the modeling of rapid chemical reactions. These systems read:

∂ t U + ∂ x F (U )= 1 ε S(f )
and are related to the so called relaxation methods (see [START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF][START_REF] Pareschi | Central differencing based numerical schemes for hyperbolic conservation laws with relaxation terms[END_REF]) that consists in replacing a nonlinear problem of the above form by a relaxed linear system of the following form:

∂ t U + ∂ x V =0,∂ t V + a∂ x U = -1 ε (V -F (U ))
where ε is called relaxation rate and a is an arbitrary velocity (with -√ a<F ′ (U ) < √ a). We remark that there is another area in the hyperbolic systems field that is related to such asymptotic namely the kinetic schemes: in this direction, one replaces the hyperbolic system of interest (e.g. the Euler equation) by a kinetic formulation in the hydro-dynamical limit (e.g. a B.G.K type equation), see [START_REF] Lions | Existence of entropy solutions for the hyperbolic systems of isentropic gas dynamics in Eulerian and Lagrangian coordinates[END_REF][START_REF] Perthame | An introduction to kinetic schemes for gas dynamics. An introduction to recent developments in theory and numerics for conservation laws[END_REF]. Then, one derives a scheme on the kinetic formulation and takes its hydro-dynamical limit. These methods have stability, accuracy and efficiency advantages. Some results on convergence are also available, see [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF].

Let us also mention another recently way to design Asymptotic Preserving Schemes, the Well Balanced schemes, see [START_REF] Gosse | Ap r i o r ierror estimate for a well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF][START_REF] Gosse | A well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF][START_REF] Greenberg | A well balanced scheme for the numerical processing of source terms in hyperbolic equations[END_REF][START_REF] Leveque | Balancing source terms and flux gradients in high-resolution Godunov methods: the quasi-steady wavepropagation algorithm[END_REF], designed to capture stationary fluxes for hyperbolic systems with source terms. In [START_REF] Gosse | Ap r i o r ierror estimate for a well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF][START_REF] Gosse | A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Gosse | A well-balanced scheme designed for inhomogeneous scalar conservation laws[END_REF] such a scheme is used for Goldstein-Taylor type models (two characteristic speeds). A Chapman-Enskog expansion shows that the asymptotic limit of the scheme is the good one for Goldstein-Taylor type models, but on systems that have more than two characteristic speeds the asymptotic seems hard to generalize. There is a large number of applications based on transport equations having a diffusive asymptotic. Let us mention, for instance, neutron transport, radiative transfer in the "optically thick limit" (see [START_REF] Adams | Subcell balance methods for radiative transfer on arbitrary grids[END_REF][START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes[END_REF][START_REF] Larsen | Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes[END_REF][START_REF] Miller | Finite differences versus finite elements in slab geometry, even-parity transport theory[END_REF][START_REF] Morel | A linear-discontinuous spatial differencing scheme for Sn radiative calculations[END_REF] and the references therein) and semiconductor modeling (see [START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Klar | An asymptotic-induced scheme for non stationary transport equations in the diffusive limit[END_REF][START_REF] Markowich | Semiconductor equations[END_REF]). However, many of the previous works deal with the so-called telegrapher equation, or equivalently Goldstein-Taylor equation which is a kinetic equation where the distribution function is localized on two opposite velocities (see [START_REF] Caflisch | Uniformly accurate schemes for hyperbolic systems with relaxation[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | Numerical integrations of systems of conservation laws of mixed type[END_REF][START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF]). On the other hand, in this paper, we shall not separate the particle density function with respect to the sign of the velocity, usually called parity method or even-odd decomposition. We refer more precisely to the introduction of Section 4.4 for a detailed explanation of this fact and a comparison with previous works. Let us finally mention that the space discretization scheme we propose could be also apply to the Boltzmann-Lorentz operator with an arbitrary cross-section.

This work is divided as follows. In Section 2, we briefly recall the Hilbert expansion method at the continuous level in order to derive the diffusion model. Then, we consider successively the discretization with respect to the velocity angle θ, to the space variable x, and finally to the time t.

In Section 3, following Jin-Levermore (see [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF][START_REF] Jin | Fully-discrete numerical transfer in diffusive regimes[END_REF], for what regards the isotropic Boltzmann-Lorentz collision operator), we consider the problem discretized only with respect to the velocity variable. We look for a choice of velocity discretization points such that the diffusion coefficient tends towards the value obtained in the continuous case. Our goal is also to use a small number (denoted by N θ ) of discretization points with respect to velocity angle, in order to avoid too much expensive computations. We prove that if we consider a uniform grid, we obtain the right diffusion coefficient only in the limit N θ →∞(with an accuracy of order 2). In order to use a kinetic description with a small number of discretization points and to preserve the right asymptotic when ε → 0, it is thus necessary to consider a non-uniform (but symmetric) grid. In particular, we treat the case of a grid with 4 or 8 discretization points.

In Section 4, we consider the space discretization (the velocity discretization is assumed to be known). More precisely, we study the discretization of the transport term v∂ x f . The boundary layer problem is not taken into account, and will be studied in futures works. The first approach consists of discretizing this term by means of an upwind scheme, but this method leads to a parasite solvability condition, i.e. we get an infinite diffusion coefficient. On the contrary, the centered scheme converges to a discretization of the Laplacian (with respect to x), but the discretization acts on a double mesh, the even points of the mesh being decoupled from the odd points, yielding to some spurious modes (which gives numerical oscillations in the computations). We then consider a θ-scheme (ε part for the upwind scheme and (1 -ε) for the centered one): the first part effectively re-couples the discretization points, but it introduces an error of order ∆x in the value of the diffusion coefficient. Finally, we take into account the modified Jin-Levermore scheme which is obtained following the strategy proposed by Jin and Levermore (see [START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | The relaxation schemes for systems of conservation laws in arbitrary space dimensions[END_REF]). It consists in writing a finite volume type scheme and in computing the fluxes at the interfaces using the leading order term (or equivalently, the steady state equation) in the upwind scheme for the half mesh. We can also give an interpretation of this scheme as a finite discontinuous elements P 1 scheme (see [START_REF] Samba | Limite asymptotique d'un schéma d'éléments finis linéaires discontinus lumpés en régime diffusion[END_REF]).

Section 5 deals with the time discretization. First, we remark that the usual methods of splitting transport and collisions are not suitable in the limit ε → 0. In fact, both the collision and the transport parts yield to projecting onto the constant states (with respect to velocity or position), and in only two time iterations one may get a constant function. Moreover, to avoid the CFL condition, which would lead to a very small time step when ε → 0, we must use implicit schemes. Numerical results are given in Section 6. Some final results and comments are then presented as a conclusion in Section 7.

The diffusion limit at the continuous level

Let us first recall the derivation of the diffusion equation from the kinetic one by means of the Hilbert expansion method. As explained in the introduction, we consider the following kinetic equation of unknown f = f (x, θ, t) in one space variable x ∈ R, the angular velocity being θ ∈ (0, 2π) with periodic condition on θ, at the diffusion scale (t/ε 2 ,x/ε):

ε∂ t f +cosθ∂ x f = 1 ε ∂ 2 f ∂θ 2 , (2.1)
where the collision term is the Fokker-Planck-Lorentz operator and ε ≪ 1 denotes the Knudsen number. We now use a classical Hilbert method by expanding f in terms of ε:

f = f 0 + εf 1 + ε 2 f 2 + ...,
and identifying terms of equal powers in (2.1). We successively obtain:

∂ 2 f 0 ∂θ 2 =0 , (2.2) cos θ∂ x f 0 = ∂ 2 f 1 ∂θ 2 , (2.3)
∂ t f 0 +c o sθ∂ x f 1 = ∂ 2 f 2 ∂θ 2 • (2.4)
From the first equation (2.2) and periodicity, we deduce that f 0 is independent of the velocity angle θ.T h e solvability condition for the second equation writes:

2π 0 cos θ∂ x f 0 dθ =0, (2.5)
and this condition is actually satisfied, because f 0 does not depend on θ and 2π 0 cos θ dθ =0 . Moreover, we have (up to the addition of a function which only depends on the space variable):

f 1 = -cos θ∂ x f 0 . (2.6)
Finally, the solvability condition for equation (2.4) writes: 2π∂ t f 0 + 2π 0 cos θ∂ x f 1 dθ = 0, which, on account of (2.6) and the fact that 2π 0 cos 2 θdθ = π gives:

∂ t f 0 - 1 2 ∂ 2 xx f 0 =0 . (2.7)
This is the limit diffusion equation we obtain at the continuous level with the diffusion coefficient ν =1/2.

The velocity discretization

In this part, we consider the discretization with respect to the velocity variable, the time and space variables remaining continuous. This is the so-called discrete ordinate method. We shall perform the same analysis as in the continuous case (see Sect. 2). Our goal is to choose a small number of discretization points such that we obtain a discrete version of the diffusion equation (2.7) with still the right diffusion coefficient, i.e. ν =1 /2. Let us first consider a uniform grid.

Uniform grid

Let us define the uniform discretization of S 1 as the angle θ j = j∆θ with ∆θ =2 π/N θ and j =0 ...N θ -1, where the indices j are denoted modulo N θ : i.e. j = N θ corresponds to the same velocity as j =0 . T h e discretized version of (2.1) is a finite system of transport equations coupled by the collision term:

ε∂ t f j +cosθ j ∂ x f j = 1 ε (Lf ) j (3.1)
where f j (x, t)=f (x, θ j ,t) is the value of the discretized distribution function f at time t, position x and velocity (cos θ j , sin θ j ). The operator L (with entries L ij ) in the right hand side of equation (3.1) stands for the discretized version of the Laplacian operator with respect to the variable θ. Using standard finite difference techniques and taking into account periodic boundary conditions, we obtain:

L ij = 1 ∆θ 2 -2i f i = j 1i f |i -j|≡1[N ].
Let us recall some useful spectral properties of the matrix L:

• The kernel of the matrix L is generated by the vector 1 =(1, 1, 1, ..., 1).

• The vector (cos θ j ) j=0...N θ-1 is an eigenvector of matrix L; the associated eigenvalue is -λ where λ is the positive number given by:

λ =2 1 -cos(∆θ) (∆θ) 2 • (3.2)
Performing the Hilbert expansion on the discretized problem (3.1), i.e.: f j = f 0 j + εf 1 j + ε 2 f 2 j + ... and identifying terms of equal powers of ε,w eg e t :

(Lf 0 ) j =0 , (3.3) cos θ j ∂ x f 0 =( Lf 1 ) j , (3.4)
∂ t f 0 +cosθ j ∂ x f 1 =( Lf 2 ) j . (3.5)
At the order ε -1 , we obtain that f 0 j does not depend on θ j , since the kernel of L consists of constant vectors; we simply set: f 0 j = f 0 . Now, the solvability condition for equation (3.4) reads:

N θ -1 j=0 cos θ j ∂ x f 0 = 0, which, because f 0 is independent of θ j ,l e a d st o : N θ -1 j=0 cos θ j =0. (3.6)
We can interpret relation (3.6) as a symmetry condition. This condition is satisfied by our choice of discretization points (with N θ even), but it does not hold if we translate the discretization points into θ 0 +2πj/N θ with θ 0 =0 for example. We look for a solution f 1 of (3.4) orthogonal to the kernel of L (Ker(L) ⊥ = {f : <f,1 >=0 }), i.e. such that:

N θ -1 j=0 f 1 j =0. (3.7)
Then, as f 0 is independent of θ j and (θ j ) satisfies (3.6), the unique solution of (3.4) reads:

f 1 j = - cos θ j λ ∂ x f 0 . (3.8)
Finally, considering (3.8), equation (3.5) becomes:

∂ t f 0 - cos 2 θj λ ∂ 2 xx f 0 =( Lf 2 ) j .
The solvability condition for this equation ( j (Lf 2 ) j =0)writes:

∂ t f 0 - 1 2λ ∂ xx f 0 =0, (3.9)
because we have: (1/N θ ) j cos 2 θ j =1 /2. We then get a diffusion equation for f 0 , but with a diffusion coefficient ν = 1 2λ : in fact, this coefficient is never equal to the expected value 1/2 (since λ = 1, see (3.2)), but it converges towards it when N θ tends to ∞ (with an accuracy of order (∆θ) 2 ). Nevertheless, if we want to deal with small number of points in θ, for lowering the computational cost, we have to investigate non-uniform grids. When using a large number of points, one can use an uniform discretization.

Non-uniform grid

First, let us briefly recall the discretization of the second derivative with respect to the angle using nonuniform grids, in order to introduce the notations. Given a discretization of S 1 , θ j for j =0 ..

.N θ -1, let us denote ∆ + j =( θ j+1 -θ j )a n d∆ - j =(θ j -θ j-1
). Applying a Taylor expansion with respect to the point θ j we obtain:

f (θ j +∆ + j )=f j +∆ + j f ′ j + ∆ + j 2 2 f ′′ j + o(∆ + j 2 ), f (θ j -∆ - j )=f j -∆ - j f ′ j + ∆ - j 2 2 f ′′ j + o(∆ - j 2 ).
By combining these Taylor expansions, we get the following approximation of the second derivative at θ j :

f ′′ j ≃ 2 f (θ j +∆ + j ) ∆ + j (∆ + j +∆ - j ) +2 f (θ j -∆ - j ) ∆ - j (∆ + j +∆ - j ) -2 f j ∆ + j ∆ - j • (3.10)
We shall discretize the Laplace operator using this formula when N θ = 4 or 8. We remark that in the general case the vector with components cos θ j is no more an eigenvector of L, as it is easily seen replacing f j by cos θ j in (3.10).

Let us perform the same discretized Hilbert expansion as presented in the uniform grid case (Sect. 3.1). At first order, f 0 j = f 0 is independent on θ j . Let us define Y ∈ R N θ as the unique solution of (LY ) j =cosθ j , for all j ∈{0, ..., N θ -1}, (3.11) such that

N θ -1 j=0 Y j =0.
Then, the unique solution f 1 of equation (3.4) satisfying (3.7) is given by:

f 1 = Y∂ x f 0 .
Equation (3.9) now reads:

∂ t f 0 +cosθ j Y j ∂ 2 xx f 0 =( Lf 2
) j , and the solvability condition for this equation of unknown f 2 still gives a diffusion equation for f 0 where the diffusion coefficient ν is now given by:

ν = - 1 N θ N θ -1 j=0 Y j cos θ j = - 1 N θ LY • Y. (3.12)
We shall now show that for N θ =4orN θ = 8, one can construct symmetric sets of discretization points that give the good diffusion coefficient. This can be related to the analysis presented in [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF][START_REF] Jin | Fully-discrete numerical transfer in diffusive regimes[END_REF] for the Boltzmann-Lorentz operator.

points of discretization

Let us first consider the case N θ =4,i.e. 4 points of discretization in velocity: θ 0 ,θ 1 ,θ 2 ,θ 3 . We still consider a symmetric discretization, i.e. such that (3.6) holds. More precisely, we choose θ 0 , θ 1 = π -θ 0 , θ 2 = π + θ 0 , θ 3 =2π -θ 0 i.e. cos θ 0 = -cos θ 1 = -cos θ 2 =cosθ 3 (see Fig. 1). The goal is to determine θ 0 in such a way that the associated diffusion coefficient ν is exactly equal to 1/2.

The discrete Laplacian matrix corresponding to this nonuniform grid reads:

L =2            -1 ∆ 0 ∆ 3 1 ∆ 0 (∆ 0 +∆ 3 ) 0 1 ∆ 3 (∆ 0 +∆ 3 ) 1 ∆ 0 (∆ 0 +∆ 1 ) -1 ∆ 0 ∆ 1 1 ∆ 1 (∆ 0 +∆ 1 ) 0 0 1 ∆ 1 (∆ 1 +∆ 2 ) -1 ∆ 1 ∆ 2 1 ∆ 2 (∆ 1 +∆ 2 ) 1 ∆ 3 (∆ 2 +∆ 3 ) 0 1 ∆ 2 (∆ 2 +∆ 3 ) -1 ∆ 2 ∆ 3           
w h e r ew ed e n o t eb y∆

0 = θ 1 -θ 0 ,∆ 1 = θ 2 -θ 1 ,∆ 2 = θ 3 -θ 2 ,∆ 3 = θ 0 -θ 3 the discretization steps.
Proposition 3.1. Given N θ =4 , there exists an unique value of θ 0 such that we have the right diffusion coefficient , i.e. ν =1/2.

Proof. We have to find

Y =(Y 0 ,Y 1 ,Y 2 ,Y 3 ) ∈ R 4 such that: LY =(cosθ)wit h 3 j=0 Y j = 0
, where the notation (cos θ) denotes the vector with components (α, -α, -α, α), with α =cosθ 0 . After some easy computations, we find

Y = -B(cos θ), with B = π(π -2θ)/4,
and the diffusion coefficient is given by ν(θ 0 )=- 1 4 Y • (cos θ)=Bα 2 . In order to find the right diffusion coefficient, we now have to look for a value θ 0 ∈ [0,π/2] such that ν(θ 0 )=1/2, where

ν(θ)= π 4 (π -2θ)cos 2 θ.
We have ν(0) = π 2 /4 > 1/2, ν(π/2) = 0 < 1/2a n dν ′ (θ 0 ) < 0. Thus, there exists a unique θ 0 ∈]0,π/2[ such that ν(θ 0 )=1/2.

Surprisingly, for the non uniform (but symmetric) discrete point, the vector with components cos θ j is an eigenvector, associated with the eigenvalue -B -1 , for the matrix L. An approximated value of θ 0 ,n a m e l y θ 0 =0 .8462 has been computed using Maple . This value differs from π/4, which is the case of uniform grid, but it is close to it.

points of discretization

We now consider the case N θ =8,i.e. 8 points of discretization in velocity. Let us denote by θ j with j =0,... ,7 a general 8-points discretization of S 1 .

We choose a symmetric discretization, i.e. such that condition (3.6) holds, as follows: θ 0 and θ 1 belonging to [0,π/2] (the others angles being induced by this choice (see Fig. 2)). We now have two angles (θ 0 and θ 1 ) that we shall determine in such a way that first the diffusion coefficient is the right one, and secondly that Y is an eigenvector (this second condition is arbitrary; we add it in order to determine the two angles in a unique way and by analogy with previous cases). The discretization step, which we recall is not fixed, is given by: We also impose that θ 1 = mθ 0 for 0 <θ 0 <π/2a n d1<m< π 2θ0 . With the same method applied in the case of 4 velocities, we can derive the discretized version of the Laplacian L:

∆ 0 = θ 1 -θ 0 , ∆ 1 = θ 2 -θ 1 = π -2θ 2 , ∆ 2 = θ 3 -θ 2 = θ 1 -θ 0 , ∆ 3 = θ 4 -θ 3 =2θ 0 , ∆ 4 = θ 5 -θ 4 = θ 1 -θ 0 , ∆ 5 = θ 6 -θ 5 = π -2θ 1 , ∆ 6 = θ 7 -θ 6 = θ 1 -θ 0 , ∆ 7 = θ 0 -θ 7 =2θ 0 .
L =2                  -1 ∆ 0 ∆ 7 1 ∆ 0 (∆ 0 +∆ 7 ) 0 ... 0 1 ∆ 7 (∆ 0 +∆ 7 ) 1 ∆ 0 (∆ 0 +∆ 1 ) -1 ∆ 0 ∆ 1 1 ∆ 1 (∆ 0 +∆ 1 ) 0 . . . 0 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 ... 0 1 ∆ 5 (∆ 5 +∆ 6 ) -1 ∆ 5 ∆ 6 1 ∆ 6 (∆ 5 +∆ 6 ) 1 ∆ 7 (∆ 6 +∆ 7 ) 0 ... 0 1 ∆ 6 (∆ 6 +∆ 7 ) -1 ∆ 6 ∆ 7                  =2             BB 1 00000 B 2 A 1 AA 2 00000 0 A 2 AA 1 0000 00 B 1 BB 2 000 000 B 2 BB 1 00 0000 A 1 AA 2 0 00000 A 2 AA 1 B 2 00000 B 1 B            
w h e r ew eh a v es e t

A = -1/αγ, A 1 =1/γ(γ + α),A 2 =1/α(γ + α), B = -1/βγ, B 1 =1/γ(γ + β),B 2 =1/β(γ + β),
with: α = π -2mθ 0 , β =2θ 0 and γ = θ 0 (m -1). We still look for a vector Y ∈ R 8 solution of (3.11) and such that it is an eigenvector for the matrix L, i.e.:

LY =(cosθ)=-λY. (3.13)
We remark that now the vector (cos θ)i sg i v e nb y( c o sθ)=(α, β, -β, -α, -α, -β, β, α), where α =cosθ 0 and β =cosθ 1 =cos(mθ 0 ). From system (3.13) we get the following system:

2B 1 (β -α)=-λα, (3.14)
-2A 1 (-α + β) -4βA 2 = -λβ.
We recall that α, β and γ depend on (m, θ 0 )s ot h a tA, A 1 ,A 2 ,B,B 1 and B 2 also depend (m, θ 0 ). Thus, equation (3.14) yields to

(π -2mθ)(π -mθ -θ)(cos(mθ) -cos(θ)) cos(mθ)=(mθ + θ)cos(θ)[(π -2mθ)cos(θ) -(π -2θ)cos(mθ)]. (3.15)
It remains to determine θ 0 and m such that the diffusion coefficient ν is equal to (1/2) i.e. such that:

-2(cos(mθ 0 ) -cos θ 0 )) = (m 2 -1)θ 2 0 (cos 2 θ 0 +cos 2 (mθ 0 )) cos θ 0 . ( 3

The space discretization

We use a finite difference discretization (see [START_REF] Godlewski | Numerical approximations of hyperbolic systems of conservation laws[END_REF]) in order to approximate the transport problem (2.1). We recall that the aim is to derive a discrete model which at the limit ε → 0 gives a good approximation of the diffusion equation (2.7). The velocity discretization is done using the finite difference scheme described in Section 3 (in particular, we shall set N θ = 4 in the simulations), and we recall that L is the discrete operator approximating the Laplacian ∂ 2 ∂θ 2 . We now describe different schemes for the convective part and study the asymptotic properties of the global space-velocity scheme. In all the sequel, we denote by f i,j an approximation of f (x i ,θ j ,.), where x i = i∆x are the points of discretization in space variable (for simplicity we set h =∆ x), and v j =( c o s θ j , sin(θ j )), j ∈{0, ..., N θ -1} are those with respect to the velocity variable. We recall in particular that these points θ j are symmetric in the sense that they satisfy the symmetry condition (3.6), which is also the discrete analogous of the continuous solvability condition (2.5). We also recall that these points are not necessarily equi-distributed, and in particular that for N θ small (N θ = 4 for example ) we have shown that, in order to recover the good diffusion coefficient, these points are necessarily distributed in a non uniform way (we still refer to Sect. 3). Furthermore, denoting by Y the vector in R N θ such that (LY ) j =cosθ j , for all j ∈{0, ..., N θ -1} with Y j = 0, the discrete diffusion coefficient is then given by (3.12) and at least in the simple cases N θ = 4 or 8 there exists a unique choice of such discrete points θ j such that ν is the good diffusion coefficient, i.e. such that:

ν = 1 N θ j Y j cos θ j = 1 2 •
In what follows, the space-velocity discretized kinetic equation reads, for all i, j:

ε∂ t f i,j +c o sθ j (Df) i,j = 1 ε (Lf i,. ) j , (4.1)
where D denotes a finite difference operator approximating the first derivative with respect to the space variable.

The upwind scheme

The more natural way to discretize the first derivative ∂ x f is the upwind scheme D u , which is defined by:

D u f i,j =          f i,j -f i-1,j h , if cos θ j > 0, f i+1,j -f i,j h , if cos θ j < 0.
We consider the discrete equation (4.1) with D = D u and use an Hilbert expansion, for all i, j

f i,j = f 0 i,j + εf 1 i,j + ε 2 f 2 i,j + ... . (4.2)
Identifying terms of equal powers, we get, for all i, j:

(Lf 0 i,. ) j =0 , (4.3)
cos θ j (D u f 0 ) i,j =( Lf 1 i,. ) j , (4.4)
∂ t f 0 i,j +c o sθ j (D u f 1 ) i,j =( Lf 2 i,. ) j . (4.5)
Equation ( 4.3) shows that f 0 does not depend on j. Now, equation (4.4) is solvable only if the following solvability condition is fulfilled:

N θ -1 j=0 cos θ j (D u f 0 ) i,j =0 .
But, on account of condition (3.6), we deduce that this is equivalent to:

f 0 i+1 -2f 0 i + f 0 i-1 =0 ,
which is the discrete approximation of a stationary diffusion equation. This is not the good diffusion limit we expected, since it gives an infinite diffusion coefficient ν = ∞. Hence, this scheme is not asymptotic preserving, since it gives a wrong diffusion equation at the limit ε → 0.

The centered scheme

Let us now examine the centered scheme which we shall denote by D c :

D c f i,j = f i+1,j -f i-1,j 2h •
The discrete problem is then (4.1) with D = D c . Identifying again the terms of equal powers in the Hilbert expansion (4.2) we get, for all i, j:

(Lf 0 i,. ) j =0 , (4.6) cos θ j (D c f 0 .,j ) i =( Lf 1 i,. ) j , (4.7)
∂ t f 0 i,j +c o s ( θ) j (D c f 1 .,j ) i =( Lf 2 i,. ) j . (4.8)
Equation (4.6) shows that f 0 still does not depend on j. The solvability condition of equation (4.7) is then automatically satisfied, on account of the symmetry property (3.6). Moreover, we have:

f 1 i,j = Y j D c f 0 i , (4.9)
where Y is given by (3.11). The solvability condition for equation (4.8) writes:

N θ ∂ t f 0 i + N θ -1 j=0 cos θ j D c f 1 .,j i =0 ,
which gives, on account of (4.9) and of the definition of ν, the following discrete diffusion equation:

∂ t f 0 i -ν ∆ 2h f 0 i =0 , (4.10)
w h e r ew eh a v eu s e dt h en o t a t i o n :

(∆ 2h φ) i = φ i+2 -2φ i + φ i-2 (2h) 2 •
In the same way, we shall denote by ∆ h the classical "three points scheme":

(∆ h φ) i = φ i+1 -2φ i + φ i-1 h 2 •
Although the discrete equation (4.10) is a consistent approximation of the continuous diffusion equation (2.7), we remark that this scheme will generate numerical oscillations, due to spurious modes: the discrete points corresponding to an even index i do not influence those corresponding to an odd one. We have actually observed this phenomena, by doing numerical experiments with an initial data of Dirac type (see Fig. 3). Since this scheme generates numerical oscillations, we drop it.

The ε scheme D

An interesting idea however is to combine this scheme to the upwind one (since the last scheme gives a diffusion operator of type ∆ h instead of ∆ 2h ), in order to avoid this phenomena. But we have seen that the upwind scheme does not give the right diffusion equation, on account of the solvability condition obtained when identifying constant terms. So a good compromise consist in studying the following convex combination:

D ε =( 1-ε)D c + εD u ;
we shall denote it by "ε-scheme".

The discrete problem is here (4.1) with D = D ε and the Hilbert expansion method gives, for all i, j:

Lf 0 i,. j =0 , (4.11) cos θ j D c f 0 i,j = Lf 1 i,. j , (4.12
)

∂ t f 0 i,j +c o sθ j D c f 1 .,j i -D c f 0 .,j i + D u f 0 i,j = Lf 2 i,. j . (4.13)
From equation (4.11), we still have that f 0 independent of j, and the solvability condition for equation (4.12) is fulfilled, on account of (3.6). Moreover, as (4.12) coincide with (4.7), we have (4.9). The solvability condition for equation (4.13) is then:

N θ ∂ t f 0 i + N θ -1 j=0 cos θ j D c f 1 .,j i + N θ -1 j=0 cos θ j D u f 0 i,j =0 ,
which gives, on account of expression (4.9), of the definition of ν and of the condition (3.6), the following discrete diffusion equation:

∂ t f 0 i -ν ∆ 2h f 0 i -ν ′ ∆ h f 0 i =0 ,
w h e r ew eh a v es e t :

ν ′ = h 2N θ N θ -1 j=0 | cos θ j |.
This equation is still a consistent approximation of the diffusion equation (2.7) and it links the points with even index to those with an odd one. Thus, we can expect this scheme to suppress the oscillations given by the pure centered scheme. This has been effectively observed numerically, but for an initial data which support contains at least two grid points (in space variable). For a Dirac type initial data however, the oscillations still persist (see Fig. 4). This scheme thus seems to be a relatively good one, although the diffusion coefficient is not exactly the good one (up to an 0(h) error). We now examine a last scheme which seems to have all the expected properties.

The modified Jin-Levermore scheme

This new scheme is inspired by the scheme proposed by Jin and Levermore for the telegraph problem in [START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF] (which is based on the use of the steady states approximation method, see for example [START_REF] Glimm | A generalized Riemann problem for quasi one dimensional gas flows[END_REF][START_REF] Vanleer | On the relation between the upwind differencing schemes of Engquist-Osher, Godunov and Roe[END_REF], ...). In this work, Jin and Levermore studied semi-discrete numerical schemes for hyperbolic systems with stiff relaxation terms that have a long time behavior governed by reduced systems of parabolic type, but the diffusion terms are in fact corrective terms of order ε. A similar problem, but under diffusive scaling, called the Goldstein-Taylor model (see [START_REF] Goldstein | On diffusion by discontinuous movements, and on the telegraph equation[END_REF][START_REF] Taylor | Diffusion by continuous movements[END_REF]), has been numerically considered by [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF]: here, both the convective terms and the relaxation terms are stiff (like for our problem), which gives additional difficulties. One of the mean ideas consists in writing the stiff terms as source terms, then to use a splitting in time algorithm, treating the relaxation terms by an implicit scheme and the new non-stiff convective terms by classical upwind schemes (or more accurate second order schemes with slope limiters). Special care has been taken to assure that the scheme possess the corrective diffusive limit. An extension of this work to more general source terms have been recently considered, in [START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF] for a bi-dimensional Boltzmann-Lorentz type operator, and for the linear Boltzmann equation in [START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF]. Both works are concerned with the actual diffusive regime and are both based on a splitting of the distribution into its even and odd part, giving then a system which is treated with similar discretization than the Goldstein-Taylor model previously described. A new difficulty arises in the Boltzmann case to solve the collision step implicitly: one needs in fact to invert an integral operator, which is not easy to do in an efficient way. A velocity discretization using Hermite polynomials is then performed.

It is in fact possible to adapt this method here to our problem, by first splitting f into an even and an odd part and then using the discrete algorithm proposed in [START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF]. Moreover, let us point out that the implicit treatment of the collision part is here far simpler for our model than for the Boltzmann case considered in [START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF]. For the continuous in time scheme, we obtain the good diffusion limit with the right coefficient of diffusion. For the full discrete problem, we obtain the right diffusion coefficient, up to an O(h) error, like for the ε scheme proposed in Section 4.3.

It is one of the reasons why we decided to construct a new scheme, still based on Jin-Levermore's approach, but where it would be possible to directly discretize f (without needing to split it into and even and an odd part) which would be far simpler. Moreover, the splitting method proposed in [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF] seems interesting in so far the discretization with respect to the space variable is "unconnected" from the velocity one, but as the stiff convective terms appear in the right hand side of the relaxation part, it seems that both meshes are coupled, giving this scheme less attractive. Finally, we point out that the scheme proposed below can be also applied to the Boltzmann-Lorentz operator with an arbitrary cross-section.

The idea for the construction of our new scheme is based on an evaluation of the fluxes at the interface between the two cells [x i -h/2,x i + h/2] and [x i+1 -h/2,x i+1 + h/2]. The scheme writes (4.1), with D = D 1/2 defined by:

(D 1/2 f ) i,j = f i+1/2,j -f i-1/2,j h •
The idea for the evaluation of the interface values f i+1/2,j (or f i-1/2,j ) is first to write (up to an O(h 2 ) error) the Taylor expansion of f i+1/2,j :

f i,j ≃ f i+1/2,j - h 2 (∂ x f .,j ) i+1/2 , if cos θ j > 0, f i+1,j ≃ f i+1/2,j + h 2 (∂ x f .,j ) i+1/2 , if cos θ j < 0,
which corresponds to a classical upwind scheme, according to the sign of the velocity. Now, in order to compute the space gradients, we consider the following continuous equation, with respect to the space variable only, i.e.

ε∂ t f +c o sθ∂ x f = 1 ε Lf,
in which we neglect the lowest order terms, with respect to ε, i.e. the O(ε) term. This gives the following system (for simplicity, we suppose that, from now on, we never have cos θ j =0):

f i,j = f i+1/2,j - h 2ε cos θ j (Lf i+1/2,. ) j , if cos θ j > 0, (4.14) f i+1,j = f i+1/2,j + h 2ε cos θ j (Lf i+1/2,. ) j , if cos θ j < 0. (4.15)
In other words, the interface values satisfy:

Id - h 2ε| cos θ j | L f i+1/2,. j = f i,j , if cos θ j > 0, f i+1,j , if cos θ j < 0.
We now use the Hilbert method, expanding f i,j according to (4.2) and identifying terms of equal powers in the kinetic equation (4.1); we get, for all i, j:

Lf 0 i,. j =0 , (4.16
)

cos θ j D 1/2 f 0 i,j = Lf 1 i,. j , ( 4
.17)

∂ t f 0 i,j +c o sθ j D 1/2 f 1 i,j = Lf 2 i,. j . (4.18)
The first equation shows, as usual, that f 0 i,j does not depend on j. Now, in order to solve the other equations, we also need to expand f i+1/2,j in terms of ε.W es e tf o ra l li, j

f i+1/2,j = f 0 i+1/2,j + εf 1 i+1/2,j + ε 2 f 2 i+1/2,j + ...
and identify terms of equal powers in equations (4.14, 4.15), this gives:

Lf 0 i+1/2,• j =0 , (4.19
)

Lf 1 i+1/2,. j = - f 0 i -f 0 i+1/2,j h/2 cos θ j , if cos θ j > 0, (4.20
)

Lf 1 i+1/2,. j = f 0 i+1 -f 0 i+1/2,j h/2 cos θ j , if cos θ j < 0. (4.21)
We first deduce from (4.19) that f 0 i+1/2,j does not depend on j; we simply denote it by f 0 i+1/2 . To compute f 1 i+1/2,. from equations (4.20, 4.21), there appears a necessary condition of solvability which writes:

-   j,cos θj>0 cos θ j   f 0 i -f 0 i+1/2 h 2 +   j,cos θj <0 cos θ j   f 0 i+1 -f 0 i+1/2 h 2 =0 , i.e. f 0 i -f 0 i+1/2 + f 0 i+1 -f 0 i+1/2 =
0, on account of (3.6). We deduce that f 0 i+1/2 is the mean value of f 0 at points x i and x i+1 , i.e. if we have

f 0 i+1/2,j = f 0 i+1/2 = f 0 i + f 0 i+1 2 •
Injecting this expression in (4.20, 4.21), we get, for all j such that cos θ j =0:

Lf 1 i+1/2,. j = f 0 i+1 -f 0 i h cos θ j ,
which gives:

f 1 i+1/2,j = Y j f 0 i+1 -f 0 i h • (4.22)
We now turn back to equations (4.17, 4.18). As (D 1/2 f 0 ) i,j is independent of j, the solvability condition of equation (4.17) is satisfied. On account of expression (4.22), we have (D 1/2 f 1 ) i,j = Y j (∆ h f 0 ) i , so that the solvability condition of equation (4.18) simply writes:

∂ t f 0 i -ν(∆ h f 0 ) i =0 ,
which is exactly the discrete diffusion equation we expected, with the good coefficient of diffusion ν: this scheme seems to be the best one, which has been confirmed by numerical tests (see Figs. 5 and6).

A finite volume interpretation of the new scheme

This new scheme can be interpreted as a generalization of the P 1 discontinuous finite elements scheme proposed by G. Samba (see [START_REF] Samba | Limite asymptotique d'un schéma d'éléments finis linéaires discontinus lumpés en régime diffusion[END_REF]) for a stationary kinetic equation of Boltzmann-Lorentz type (see also [START_REF] Larsen | The asymptotic diffusion limit of discretized transport problems[END_REF][START_REF] Morel | A linear-discontinuous spatial differencing scheme for Sn radiative calculations[END_REF]). But it can also be interpreted in terms of a classical finite volume scheme expressed on a refined mesh of size h/2. We now detail this interpretation.

The aim is to explain the computation (4.14, 4.15) of the numerical flux. We first set i

′ = i -1/2a n d introduce the new discretization points x i ′ = i ′ h and x i ′ +k/4 =(i ′ + k/4)h,f o rk ∈ Z.W es e t : f i+1/2,j = 1 2 [ f i ′ +3/4,j + f i ′ +5/4,j ],
where f i ′ +3/4,j and f i ′ +5/4,j solve the following equations:

ε∂ t f i ′ +3/4,j +c o sθ j F i ′ +1,j -F i ′ +1/2,j h/2 = 1 ε (Lf i ′ +3/4,. ) j , (4.23)
ε∂ t f i ′ +5/4,j +c o sθ j F i ′ +3/2,j -F i ′ +1,j h/2 = 1 ε (Lf i ′ +5/4,. ) j , (4.24)
with:

F i ′ +1/2,j = 1 2 [ f i ′ +1/4,j + f i ′ +3/4,j ], F i ′ +1,j = f i ′ +3/4,j if cos θ j > 0, f i ′ +5/4,j if cos θ j < 0.
In the same way, we also naturally have:

ε∂ t f i ′ +1/4,j +c o sθ j F i ′ +1/2,j -F i ′ ,j h/2 = 1 ε (Lf i ′ +1/4,. ) j . (4.25)
Now equations (4.23) and (4.24) give by addition:

ε∂ t f i+1/2,j +c o sθ j F i ′ +3/2,j -F i ′ +1/2,j h = 1 ε (Lf i+1/2,. ) j ,
or equivalently:

ε∂ t f i+1/2,j +c o sθ j F i+1,j -F i,j h = 1 ε (Lf i+1/2,. ) j .
We here recover the equation of evolution (4.1) expressed at the new spatial grid point x i+1/2 . It remains to find an equation for the flux F i,j at the interface of the refined mesh. This one is obtained by summing equations (4.23-4.25); we get:

ε∂ t F i,j +c o sθ j F i ′ +1,j -F i ′ ,j h = 1 ε (LF i,. ) j . (4.26)
A simple computation first shows that:

F i ′ +1,j -F i ′ ,j h =        2 h [F i,j -f i-1/2,j ]f o r c o s θ j > 0, 2 h [f i+1/2,j -F i,j ]f o r c o s θ j < 0.
Now, neglecting the lowest order term (with respect to ε) in (4.26), we get:

2cosθ j h [F i,j -f i-1/2,j ]= 1 ε (LF i,. ) j , if cos θ j > 0, 2cosθ j h [f i+1/2,j -F i,j ]= 1 ε (LF i,. ) j , if cos θ j < 0,
or equivalently, for i ′ = i +1/2:

f i ′ ,j = F i ′ +1/2,j - h 2ε cos θ j (LF i ′ +1/2,. ) j , if cos θ j > 0, f i ′ +1,j = F i ′ +1/2,j + h 2ε cos θ j (LF i ′ +1/2,. ) j , if cos θ j < 0,
i.e. we exactly recover formulae (4.14, 4.15) at the new grid point x i ′ .

The time discretization

We now present the discretization in time. We first remark that it is of no use to apply a splitting in time method. In fact, for the first time step, the collision part

ε∂ t f = 1 ε Lf
would project on the constants with respect to the velocity angle θ, and the transport part

ε∂ t f +cosθ∂ x f =0
would project on the constant with respect to the position variable x. Thus, in two time steps we will find a constant function (and not the good diffusion equation). The fact that splitting also fails in the diffusion limit is not surprising and we also refer to the Section A.2 where a similar problem occurs for the directional splitting in the multi-dimensional case.

On the other hand, when using an explicit scheme, we get time step restrictions. Due to the diffusive scaling of the equation (1.1), the time step stability condition associated with the transport term is of the form ∆t ≤ ε∆x, (5.1) since the velocities are of modulus smaller than 1, whereas the time step condition for the collision part reads:

∆t ≤ τε 2 (∆v) 2 , (5.2)
where τ is the collision time (equal 1 in this paper). Since the aim of the studied schemes is to be used for arbitrary small values of ε, the cost of such explicit schemes will be prohibitive especially due to the collision part. Therefore, we shall always use implicit schemes for the collision part. For the transport part, we can use either fully explicit scheme (with the restriction given by (5.1)), semi-implicit scheme (for example using the ε-scheme described in Sect. 4.3) or a fully implicit scheme. More precisely, in the semi-implicit scheme, we implicit the collision and (1 -ε)D c , the first part of the transport, and treat explicitly the upwind transport part εD u . This leads us to a stability condition of the form ∆t ≤ ∆x.

From numerical point of view, the semi and fully implicit scheme requires to invert N x × N θ matrices where N x (respectively N θ ) is the number of point in the discretization with respect to x (respectively v).

Numerical results

The numerical tests are devoted to verify that the proposed discretization gives the right diffusion coefficient. We consider an initial data equal to a Dirac measure in x and uniform in θ or concentrated in θ (but this latter choice does not change the solutions when ε → 0): For such initial data, the solution f ε of problem (1.1) behaves when ε → 0 as the solution f 0 of the heat equation, which is given by:

f (x, θ, t =0)=δ x=0 or f (x, θ, t =0)=δ x=0 δ θ=0 .
f 0 (x, t)= 1 √ 4πνt exp - x 2 4νt .
In particular, the second moment of the solution is close to that of a Gaussian

f ε (t, x, θ)x 2 dx dθ ∝ 2νt.
In other words, the evolution of the second moment of f ε becomes linear in time as ε → 0 with a slope related to the diffusion coefficient of the limiting diffusive equation.

In the results presented here, we use the four velocities case presented in Section 3.2 for a given angle θ 0 . Here, and also in the following figures, ε =0 .001, the numerical scheme is implicit in time, ∆x = 40, ∆θ =4 , ∆t =0.01 and time equal to 1 (i.e. 100 iterations).

In Figure 3, we plot the particle density n(x, t), n(x, t)= f (x, θ, t)dθ, for the centered scheme. In particular, there are underlined the numerical oscillations due to the non-coupling between the odd and even meshes. Finally, the dotted line is a post-treated curve giving the mean values of n(x, t) between two successive meshes. The density n(x, t)f o rt h eε-scheme is plotted in Figure 4. We remark that the numerical oscillations in Figure 3 are smoothened by the effect of the upwind scheme. The difference between the diffusion coefficient in the ε-scheme and the centered scheme is not remarkable due to the fact that h =∆x is small. The continuous line corresponds to the "optimal" angle computed for the nonuniform grid with four points, while the dotted line corresponds to the density n(x, t) computed for the angle θ 0 = π/6.

In Figure 5, we plot the density n(x, t) for the modified Jin-Levermore scheme. The numerical oscillations have disappeared, and the diffusion coefficient is exactly the same that in the centered scheme. The continuous line corresponds to the "optimal" angle computed for the nonuniform grid with four points, while the dotted line corresponds to an angle equal to π/6. In Figure 6, we compare the two densities in the case of the centered scheme (the means value curve) and in the case of the modified Jin-Levermore scheme. We note that the Dirac mass has diffused in the same way, as it was attended from the previous analysis.

In Figure 7, we plot the curves of the diffusion coefficient ν(θ) computed from the second moment, when θ is varying. We compare the result found for the exact solution of the heat equation with the one obtained by means of the modified Jin-Levermore scheme. Both the curves intersect the value 1/2 for the optimal value of the angle θ 0 =0.864.

In Figure 8 we plot the density n(x, t) for four different values of ε at time 1. For small values of ε (0.01 and 0.00001) the curves are indistinguishable. For "large" values of ε (larger than 1) , the solution is very close to the initial data. For intermediate values (ε =0 .5), one observes the splitting of the initial delta measure into two picks: one moving to the right the other to the left, as expected in such kinetic model.

Finally, in Figure 9 we plot the diffusion coefficient curves in function of ε =0.01,... ,0.1 for three different values of the discretization angle θ 0 = π/6,π/3, 0.864. We note that for θ 0 =0.864 (which is the angle computed by Maple for the non-uniform four points discretization) the diffusion coefficient ν is close to 1/2whenε =0.01.

Conclusions

We first remark that in the ε scheme proposed in Section 4.3, it is possible to replace the upwind scheme by the modified Jin-Levermore scheme of Section 4.4, that has all the required properties for ε ≤ h. On the other hand, if ε ≥ h, then the upwind scheme seems sufficient. Thus, the ultimate choice for the discretization in space is, for example,

D =max(0, 1 -ε/h)D 1/2 + min(1,ε/h)D u . (7.1)
Let us emphasize that every method leading to the good diffusion equation (see for example [START_REF] Gosse | A well-balanced scheme using non-conservative products designed for hyperbolic systems of conservation laws with source terms[END_REF][START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Jin | Diffusive relaxation schemes for multiscale discrete-velocity kinetic equations[END_REF]), requires the inversion of a huge (but sparse) matrix, of size N x × N θ where N x is the number of discretization points in space and N θ is the number of discretization points in velocity. In fact, one needs implicit schemes in order to avoid too restrictive time step condition in terms of ε, and a splitting of space and velocity is not possible. Indeed, for all these methods in order to compute the fluxes we must solve a stationary problem (see for example

v •∇ x f 0 = L f 1 .
Such a problem is necessarily multi-dimensional, i.e. one cannot split space directions x and velocity variables v. We also remark that there is no non-trivial hydro-dynamical limit (see Sect. A.1). Moreover, the bi-dimensional case with respect to the space variable can be treated by a similar approach but the directional splitting does not give the good diffusion coefficient (see Sect. A.2). Nevertheless, the finite volume approach described in Section 4.5, could be a way to design a fully multi-dimensional asymptotic preserving scheme.

Some questions remain open:

-Take into account boundary condition (or equivalently bounded region in x): as in [START_REF] Golse | The convergence of numerical transfer schemes in diffusive regimes I: discrete-ordinate method[END_REF][START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF][START_REF] Jin | Fully-discrete numerical transfer in diffusive regimes[END_REF][START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF], the solution has a boundary layer as ε → 0 which is treated using extrapolation length. -Generalize to a non uniform grid in space and analyze the two dimensional case (in space) with unstructured meshes, or other geometries with for example spherical or axial symmetry.

-Consider inelastic collisions or force fields that couples the kinetic energy level i.e. consider a distribution function depending of the kinetic energy. -Numerical test with regions where transport dominates (ε ≫ 1) and other highly collisional (ε ≪ 1) are needed in order to validate the choice proposed in (7.1). -Study the positiveness and other properties of the solution when using Hilbert expansion method.

A. Appendix: some final remarks A.1. Hydro-dynamical limit Consider the Boltzmann kinetic equation at a hydrodynamic regime:

∂ t f +cosθ∂ x f = 1 ε L(f ).
Is it possible to make the same analysis as for the diffusion regime. We begin by the continuous level with L(f )=∂ 2 θθ f . Expanding f in powers of ε, and identifying the equal orders, we get: 0=L f 0 , (A.1)

∂ t f 0 +cosθ∂ x f 0 = L f 1 , (A.2) ∂ t f 1 +cosθ∂ x f 1 = L f 2 . (A.3)
Then solving equation (A.1) gives that f 0 is independent on the θ variable: f 0 = f 0 (x, t). The solvability condition for equation (A.2) reads:

2π∂ t f 0 + 2π 0 cos θ∂ x f 0 dθ =0,
which implies f 0 independent on t,t o o . T h u s ,f 0 = f 0 (x)a n df 1 = f 1 (x, θ)=-cos θ∂ x f 0 . Finally, the solvability condition for equation (A.3) is given by: 2π 0 -(cos θ) 2 ∂ 2 x xf 0 =0, which first gives f 0 = C ′ x + C but, thanks to the behavior at infinity, we finally get f 0 = C; we also deduce that f 1 = 0. Thus, it shows that the hydro-dynamical limit is trivial.

A.2. Bi-dimensional case

Let us consider a two dimensional case. We shall show that a directional splitting is not suitable since it leads to a wrong diffusion coefficient. We shall present the analysis on the simplest, centered case.

We denote by i the index for the x position variable, by k the index for the y position variable. For simplicity, we shall consider only the discretization with respect to the position variables x and y, keeping the velocity θ and the time t at a continuous level. Applying a centered scheme to the kinetic equation, leads to:

ε∂ t f i,k +cosθD 1 f i,k +sinθD 2 f i,k = 1 ε Lf i,k , (A.4)
where

D 1 f i,k = f i+1,k -f i-1,k 2∆x
and

D 2 f i,k = f i,k+1 -f i,k-1 2∆y •
Performing the Hilbert expansion method with

f i,k = f 0 i,k + εf 1 i,k + ε 2 f 2 i,k ,
we get at the order ε -1 :

f i,k = f 0 i,k = ρ i,k , independent on θ.
Then at the order ε 0 we obtain: cos θD 1 f 0 i,k +sinθD 2 f 0 i,k = Lf 1 i,k , (A.5) which gives:

f 1 i,k = -cos θD 1 f 0 i,k -sin θD 2 f 0 i,k .
Finally, at the order ε 1 ,w eh a v e :

∂ t f 0 i,k +cosθD 1 f 1 i,k +sinθD 2 f 1 i,k = Lf 2 i,k
and integrating with respect to θ and replacing f 1 by (A.5), this yields to the diffusion equation on a double mesh:

∂ t ρ i,k - 1 2 ρ i+2,k + ρ i-2,k -2ρ i,k (2∆x) 2 - 1 2 ρ i,k+2 + ρ i,k-2 -2ρ i,k (2∆y) 2 =0.
Still this discretization has the same problems of the one in the one-dimensional case: it does not couple the odd and even grid points of the mesh. We remark also that the same analysis of the modified Jin-Levermore scheme can be performed when using a rectangular structured mesh but it is much more complicated with a general mesh, see [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF][START_REF] Hermeline | Two coupled particle-finite volume methods using Delaunay-Voronoï meshes for the approximation of Vlasov-Poisson and Vlasov-Maxwell equations[END_REF][START_REF] Samba | Limite asymptotique d'un schéma d'éléments finis linéaires discontinus lumpés en régime diffusion[END_REF] for works in this direction. Finally, it seems that a directional splitting (alternative directions) for the modified Jin-Levermore scheme leads to a wrong diffusion coefficient. Indeed, (A.4) can be splitted in two parts:

ε∂ t f i,k = -cos θD 1 f i,k + 1 2ε Lf i,k + -sin θD 2 f i,k + 1 2ε Lf i,k .
The first bracket corresponds to the derivative with respect to x and half of the collision, the second bracket corresponds to the derivative with respect to y and half of the collision. Then, using a splitting, one solves the first part which gives when ε tends to zero, the diffusion equation in x with coefficient ν x =1/4, instead of 1/2. Thus, in order to recover the right diffusion coefficient one may not use a directional splitting.

Figure 1 .

 1 Figure 1. 4 discretization points in velocity.

Figure 2 .

 2 Figure 2. 8 discretization points in velocity.

Figure 3 .

 3 Figure 3. Centered scheme.

Figure 5 .

 5 Figure 5. Modified Jin-Levermore scheme.

Figure 8 .Figure 9 .

 89 Figure 8. Density for three values of ε.


S Zhu

M Pressigout

M Serviè Re

L Morin

G Moreau

Skyline Matching: A robust registration method between Video and GIS

Camera pose estimation from an image has long been a very active research for Video/GIS registration. Different from traditional methods which are based on direct building feature matching, this paper presents an alternative method which uses an overall Skyline feature to determine the camera's orientation. Through our skyline matching method we find out a way to transform an image/model matching problem into a curve matching problem which significantly reduced the computation complexity and make it a possible solution for real-time applications.

Introduction

With the development of digital techniques and information industry, the requirement for more detailed and updated geographic information has been growing rapidly in last two decades. Based on this demand Geographic Information System (GIS) technology is mushrooming and also widely used, especially in the fields of navigation, territory management and environment evaluation. As one result of this mutual development, new technologies like Computer Vision and Augmented Reality (AR) are integrated into GIS. At the same time 3D models are replacing traditional 2D data.

Under this context, Video/GIS registration technique is proposed as a solution for both Augmented Reality Application [START_REF] Langlotz | Sketching up the world: in situ authoring for mobile Augmented Reality[END_REF] and 3D texture update problem [START_REF] Sourimant | Toward automatic GIS-video initial registration[END_REF].

According to [START_REF] Colleu | Automatic Initialization for the Registration of GIS and Video Data[END_REF], Video/GIS registration is a fusion between videos and geo-referenced building models which are extracted from GIS. The camera pose estimation problem is really the essential challenge of Video/GIS registration issue. If both the camera location in the world coordinate system and its orientation can be determined, then it is easy to find the relationship between an image pixel 3U3D2012 and an object from the real world (such as building facades, edges, windows or some other particular features which can be used to identify such building in the image).

In the general case, a camera pose can be expressed by a projection matrix P between the 3D world coordinate system and the 2D pixel coordinate system. This projection matrix can be decomposed into three main components: [START_REF] Langlotz | Sketching up the world: in situ authoring for mobile Augmented Reality[END_REF] Where K is the camera intrinsic parameter matrix, R is the rotation matrix and t is the translation vector. In normal case, a rough translation vector which is the position of the camera can be easily acquired by GPS signal and the intrinsic parameters can also be determined by simple camera calibration technique [START_REF] Tsai | A Versatile Camera Calibration Technique for High-Accuracy 3D Maching[END_REF]. So the most difficult part of the camera pose estimation is the calculation of the rotation, which is defined by the rotation matrix R.

One of the traditional ways to solve this problem is try to calculate the projection matrix directly through several corresponding features between an image and the GIS model. However these kinds of methods usually use iterative algorithms to establish the correspondence, which is quite computationally expensive and heavily depends on the quality of features database and on features extraction in the image. The method proposed here exploits the fact that most images taken in urban environments always contain a skyline which is a border between buildings and the sky and we assume that it is quite unique for different orientation of one known position. Moreover, it can be extracted from the image much more easily and robustly than features such as buildings edges. Our method aims to present a computationally inexpensive, simple and robust method for automatic registration.

The remainder of this paper is structured as follows: first, we summarize the related previous works. Next, section 3 describes the method of skyline extraction. Section 4 explains and applies skyline matching algorithm. Results are presented in section 5 followed by the conclusions and future work. Two hypotheses are used: 1. Horizontally rectified images; 2. Calibrated camera.

Related works

Video/GIS registration technique has been growing rapidly in recent years. The essential of Video/GIS registration is aligning each pixel in the image with 3D coordinates in the GIS model. In order to achieve this registration, we need to know the exact position and orientation of the camera (the camera pose) when the picture is shot.

The most common method consists in first trying to build up a feature database of searching areas then extract specific features from the real image then searching the database to find the best correspondence. The database may either be composed of geo-referenced images like [START_REF] Zhang | Image Based Localization in Urban Environments[END_REF] or of unique features extracted from images such as in the work of Arth et al. [START_REF] Arth | Wide area localization on mobile phones[END_REF]. This kind of methods takes a lot of time and requires storing resources to establish such database and it still needs a remarkable computation power to conduct the searching and matching process. The advantage of this kind of methods is that they do not need the initial position or rough pose information.

3U3D2012

Another type of methods is focused on directly using the geometry features of the GIS model to establish the correspondence. Such a method, proposed by Reitmayr & Drummond [START_REF] Reitmayr | Going out: robust model-based tracking for outdoor augmented reality[END_REF], uses a rough textured GIS model as matching database and applies an edge-tracking system to calculate the camera pose. An improved work was done by Sourimant et al. [START_REF] Sourimant | Toward automatic GIS-video initial registration[END_REF] where a RANSAC algorithm is applied to find out the best correspondence of line features between GIS model and image. Both their works are heavily dependent on a robust initial rough pose which is provided either through differential GPS or other additional sensors such as gyroscope. The work of Bioret et al. [START_REF] Bioret | Urban Localization based on Correspondences between Street Photographs and 2D Building GIS Layer[END_REF] provides a different point of view to solve the problem by only using the geometrical information reconstructed from the image such as the angle between the facades as a query key to find the best correspondence in the GIS database. However the robustness is not perfect because ambiguity still remains in most of the cases. Direct 3D reconstruction from image has also been considered as a solution of the Video/GIS registration problem. Some works such as [START_REF] Pollefeys | Detailed Real-Time Urban 3D Reconstruction from Video[END_REF] show the possibility to completely reconstruct a 3D building within the local camera referenced system but it needs more than one camera as input source.

On account of all the limitations of previous research, we wish to develop a method which is computational inexpensive and can calculate the camera pose from a single image. Based on such requirements we need to find a feature which can represent the overall geometric character of the image and at the same time avoid direct building feature matching. Under such concern, the skyline feature becomes a perfect candidate for our method. The skyline is easier to extract and more robust for two reasons. First, there is almost no artificial object in the sky and the color of the sky is quite uniform that makes it suffers less noise during extraction. Secondly, in the general case, all occlusions happen under the location of the skyline, which makes it a feature more robust than others. Up to our knowledge, skyline features have not been employed for pose estimation from images.

However, in the robotics and automation community the characteristics of the skyline have been investigated. Ramadingam et al. [START_REF] Ramalingam | SKYLINE2GPS: Localization in Urban Canyons using Omni-Skylines[END_REF] use an upward pointing catadioptric (i.e., fisheye lens) camera and 3D urban models for precise localization without GPS using the skyline. A. Nuchter [START_REF] Nüchter | Skyline-based registration of 3D laser scans[END_REF] uses skyline features in his registration work of 3D laser scans.

Skyline extraction

The proposed method of camera pose recognition through skyline matching is under two hypotheses:

1) the image is horizontally rectified; 2) there exists one significant sky region on the image. The idea is to extract the skyline from an image with known camera position. At the same time, a panoramic (360°) skyline around a point of view (which is the camera position of the real image) is extracted from the virtual world created from the GIS model. Then, we can compare these two skylines and find a correspondence between them, so that the partial skyline from the real image does match a sub-part of the panoramic skyline from the GIS model. The first step is the extraction of these two kinds of skyline.

3U3D2012

Coordinate transformation

In order to describe the skyline feature, we need to define a coordinate system which can best represent such information. The basic idea is trying to transform the 3D coordinate of the skyline point into a 2D presentation and such presentation must be independent of the scale factor. Under such concern, we use only the Azimuth angle φ and the Elevation angle θ to describe the point P on a skyline (See Figure 1). It is indeed the so-called cylindrical coordinate system and it is a 2D coordinate system. So for any point P(x,y,z) in the Cartesian coordinate system the transformation is as follows. First, we can get the distance ρ between P' and C where P' is the projection of point P in the horizontal plane which go through the camera center C(x 0 ,y 0 ,z 0 ): √

Then we can got the Azimuth angle:

{ ( ) ( ) (2)
Then, we can calculate the elevation angle θ by using the distance between P and C: √ and their elevation differences:

(3)

Panoramic skyline

In this section we will explain how the panoramic (360 o ) skyline is generated around a known position from the GIS model. First four joint images were created (rendered) around such the same viewpoint by using OpenGL (see Figure 2).

Secondly the depth map of each image is calculated. Based on depth information of each point, the skyline points are extracted from the image together with their 3D coordinates. Finally the coordinate transformation is conducted as presented in section 3.1 and the panoramic view of the skyline is obtained (see Figure 4.B). The same method can also be used to extract the partial skyline from a synthetic image. The extraction of the skyline from a real image is much more complicated, because there is no depth information and the sky region is not as identical as in the virtual world. In theory, the skyline is characterized of a curve which: (1) is composed of edge points that separate the terrain objects and the sky, (2) locates at the upper parts of the image and extends from one side to the opposite side of the image boundary [START_REF] Lie | A robust dynamic programming algorithm to extract skyline in images for navigation[END_REF]. In order to extract the skyline we need to distinguish the sky region from other part of the image. Generally the intensity of the sky region is more uniform than in other parts of the image, which means the intensity gradient inside the sky region is small. This characteristic allows us to use a very fast and robust method, namely adaptive threshold, to find the boundary of the sky (see Figure 3.a). the panoramic skyline using the real camera internal parameters. One remark is that all non-skyline intersections are set to zero and labelled as "no-sky" in order to ignore their influence during the skyline matching process.

Fig. 1 .

 1 Fig. 1. The transformation between the Cartesian coordinate and the Cylindrical coordinate; C, the camera center, is the point of view and P is an arbitrary point in 3D space.

3U3D2012Fig. 2 .

 2 Fig.2. Create four joint images around one point of view

Fig. 3 .

 3 Fig.3. This image show the extracted sky line (green) by using adaptive threshold method Besides, a morphologic filter closing operation on original image is also used to reduce the artefacts caused by thin structures such as: telegraph poles, wires and road lights. Then, each extracted skyline point is also transformed into the same 2D coordinate system (see Figure 4.b) as

Skyline matching method

We now have both the panoramic skyline (B) from GIS and the partial skyline (A) from the real image (See Figure 5.A and 5.B for an example). Both of them are assumed to be extracted from the same camera position. In order to find their correspondence, we rely on the similarity of their shapes.

As they are already transformed into the same scale irrelevant 2D coordinate metric, they are supposed to have the same energy level. This fact simplifies the comparison process. We use a similarity factor based on the SAD (sum of absolute differences), with an outlier removal procedure, to measure the correspondence between A and B of each azimuth angle within A as follows:

Where represents partial skyline A and represents panoramic skyline B. We call equation (4) the Cross-Similarity function expressed as . For outlier removal, points with similarity below 75% are considered as outliers and discarded from the similarity evaluation. Its purpose is to make the correspondence peaks sharper and more significant. This method is under the assumption that if the recent searching point is the corresponding point of two skylines, then their energy difference (difference of elevation angle) should be less than 25%. Figure 4 shows one result of this Cross-Similarity method process with one partial skyline from a real image and one panoramic skyline generated in the same position. From this result we can find a quite significant corresponding peak which can help us to determine the camera's orientation in the horizontal plane.

Experiments and Results

In order to test the accuracy and robustness of our method, we conducted several experiments. First we tested the accuracy of the algorithm by using a set of synthetic images which is generated by 3U3D2012 using the GIS model and known projection matrices. The result is quite satisfying; our method can correctly recognize the camera orientation as long as we have more than 10% skyline appearance along the image width. Then we applied our method in several images which are taken in the city center of Rennes, France. During these experiments we not only compared the partial skyline extracted from the image with the panoramic skyline generated in the position provided by GPS.

Instead we compare it with a set of panoramic skylines located around this position in order to pick out the position with the highest similarity value, through which we not only estimated the orientation of the camera but also refined the translation vector. Figure 5 shows two results we got from our experiment, from which we can see the GIS model (blue wireframe) is well aligned with the buildings. As we can see from the result that the upper boundary of the building is better registered than the vertical boundary, which is because the algorithm is only focused on getting the best matching between skylines which is only related to the upper boundary of the building. Such characteristic also makes this method fail when the upper boundary of the building is quite different from the GIS model (as our GIS model does not contain any roof slope information it always fails when buildings have significant roof height) or such boundary is sheltered by trees or something else.

Conclusion

Based on the concern of developing an efficient and pure geometry solution for Video/GIS registration, we investigated the possibility to use skyline feature as a bridge to establish the 3U3D2012 correspondence between real images and GIS models. Through our skyline matching method we find out a way to transform an image/model matching problem into a curve matching problem which significantly reduces the computation complexity (time for one Cross-Similarity matching is around 30ms and for a complete skyline matching without optimization is around 160ms in a PC platform with Inter Core2 CPU E8400. 3.00GHz) and makes it a possible solution for real-time applications.

Its accuracy and robustness was tested on a set of both synthetic and real images. Despite its simplicity, the skyline matching method gives promising results in the case of an artificial environment like a city center. Our future work will be focused on improving its robustness in case of occlusions between building and sky such as tall trees and also improving its efficiency by using GPU calculation. On the other hand we will try to ally this method with vanishing point estimation which can be used to rectify the image into horizontal position which is the basic demand of this method.