Nicolas Crouseilles

Mohammed Lemou

An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion

An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits.

Introduction

For the description of plasma, the fundamental kinetic model is the Vlasov equation satisfied by the distribution function f (t, x, v) which depends on time t ≥ 0, on the position x ∈ R d and on their velocity v ∈ R d , with d the dimension of the problem. Even if fluid response of the plasma requires low computational efforts, the system may be far from equilibrium, and this requires the use of kinetic models.

The modelling and simulation of plasmas at the kinetic level is known to be a challenging problem due to the high dimensionality of the phase space (6 dimensions plus time for the general case), but also due to the existence of multiscale phenomena that occur. Among them, the Debye length or the mean free path which can be very small compared to the size of the device ; in the same way, the relaxation time may also introduce an additional time scale. From a numerical point of view, a classical explicit scheme has to solve these micro-scales in order to remain stable. This makes these schemes extremely time consuming.

Then reduced models are usually proposed to describe specific regimes in specific regions of the simulated device. However, such domain decomposition approaches (see [START_REF] Arnold | Low and high-field scaling limits for the Vlasov-and Wigner-Poisson-Fokker-Planck system[END_REF][START_REF] Bourgat | Coupling Boltzmann and Euler equations without overlapping, Domain decomposition and Engineering[END_REF][START_REF] Klar | Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductors equations[END_REF][START_REF] Crouseilles | A hybrid kinetic/fluid model for solving the gas dynamics Boltzmann-BGK equation[END_REF][START_REF] Crouseilles | A hybrid kinetic-fluid model for solving the Vlasov-BGK equations[END_REF]) need heavy development to for example connect the different models at the interfaces and to distinguish the fluid region from the kinetic ones at different times. The kinetic and fluid models are solved simultaneously on the different subdomains of the computational domain. Obviously, the high dimensional cases become technical and difficult since the interfaces become curves or surfaces. Then, it seems interesting to develop alternative approaches which can switch automatically from one regime to another.

In this paper, we are concerned with the development of numerical scheme for solving the Vlasov equation with collisions that are uniformly stable along the transition from kinetic to macroscopic regime. More precisely, two asymptotics are considered here: the diffusion and the high-field limits. For the first limit, the initial model is the Vlasov-Poisson-BGK model (for which the diffusion limit has been studied mathematically in [START_REF] Poupaud | Diffusion approximation of the linear semiconductor Boltzmann equation[END_REF] and numerically in [START_REF] Klar | A Numerical Method for Kinetic Semiconductor Equations in the Drift Diffusion Limit[END_REF]) whereas for the second limit, a more suitable model is the Vlasov-Poisson-Fokker-Planck equation (see [START_REF] Nieto | High-field limit for the Vlasov-Poisson-Fokker-Planck system[END_REF][START_REF] Bostan | Electric high-electric field limit for the Vlasov-Maxwell-Fokker-Planck system[END_REF][START_REF] Bonilla | High-field limit of the Vlasov-Poisson-Fokker-Planck system for different perturbation methods[END_REF] for a mathematical study). It is well known that standard discretizations suffer from restrictive condition on the time step ∆t (∆t = O(ε 2 ) or O(ε), where ε is the small parameter) which prevents from a numerical asymptotic study.

To construct such numerical schemes, we propose an extension of the micro-macro decomposition based strategy introduced in [START_REF] Liu | Boltzmann equation: micromacro decompositions and positivity of shock profiles[END_REF][START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Jin | A micro-macro decomposition based on asymptotic-preserving scheme for the multispecies Boltzmann equation[END_REF] to the framework of collisional Vlasov-Poisson kinetic equation. The micro-macro decomposition is simply based on the fact that the kinetic unknown f is equal to the sum of an equilibrium and a rest g (see Figure 1) as for the Chapman-Enskog expansion. The associated model is then composed of a kinetic equation satisfied by g and a macroscopic equation satisfied by the moments of f (which are equal to those of its associated equilibrium). The model is solved with a suitable semi-implicit method which overcomes the difficulty induced by the stiff terms when the small parameter becomes close to zero. Obviously, this is of great interest since standard explicit schemes lead to very expensive numerical simulations (the time step has to be of the order of the small parameter). One other interesting property for this approach is the following: a general class of numerical schemes can be constructed from this micro-macro decomposition that enjoys the Asymptotic Preserving (AP) property (see [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF]): if we fix the numerical parameters, and let the small parameter going to zero then the micro-macro decomposition provides a consistent numerical scheme for the asymptotic limit.

Previous works have been performed in this direction but using other strategies or applied to other models (radiative transfer, semi-conductors or Boltzmann). We can cite for the diffusion regime [START_REF] Jin | The discrete-ordinate method in diffusive regimes[END_REF][START_REF] Jin | Numerical schemes for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Jin | A micro-macro decomposition based on asymptotic-preserving scheme for the multispecies Boltzmann equation[END_REF][START_REF] Klar | Asymptotic-induced domain decomposition methods for kinetic and drift diffusion semiconductors equations[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF][START_REF] Gosse | Asymptotic-preserving and well-balanced schemes for radiative transfer and the Rosseland approximation[END_REF], especially [START_REF] Klar | A Numerical Method for Kinetic Semiconductor Equations in the Drift Diffusion Limit[END_REF] in which the diffusion limit of charged particles is investigated. Fluid limits have also been studied in [START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Filbet | A class of asymptotic preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Degond | Macroscopic fluid models with localized kinetic upscaling effects[END_REF] and we can mention works devoted to the quasineutral limit [START_REF] Degond | Asymptotic-Preserving Particle-In-Cell method for the Vlasov-Poisson system near quasineutrality[END_REF][START_REF] Belaouar | An asymptotically stable semi-Lagrangian scheme in the quasi-neutral limit[END_REF]. The main purpose is to develop a model (equivalent to the original one) the discretization of which is able to describe the different regimes that can occur in a realistic simulation. Standard asymptotic preserving schemes usually capture the limit ε = 0 but the micro-macro approach enables to capture the first order corrective terms in ε (see [START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF]). This property can be explained by the fact that the micro-macro decomposition mimics in a certain way the Chapman-Enskog method. Let us mention that some strategies are based on a splitting procedure but this is not the case here. Moreover, the present approach enables an implicit treatment of the space-diffusion term which appears in the asymptotic model (in both diffusion and high-field asymptotics), as suggested in [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF]. This improvement of the micro-macro approach has the following important property : it makes the numerical scheme free from the usual restrictive CFL condition ∆t = O(∆x 2 ) at all regimes where ∆x denotes the spatial step. Finally, as in [START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | A Numerical Method for Kinetic Semiconductor Equations in the Drift Diffusion Limit[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF], we are interested in the study of boundary conditions. This is done following [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. In this work, the starting model is the collisional Vlasov-Poisson equation coupling a transport equation on the distribution function f (t, x, v) ≥ 0 in the phase space, with the self consistent electric field E(t, x). The spatial variable x belongs to the interval [0, L], L > 0, whereas the velocity v belongs to R and the time t to R + . More precisely, the system reads

∂f ∂t + v ε α ∂ x f + E ε ∂ v f = 1 ε 1+α Q(f ), (1.1)
∂ x E(t, x) = ρ(t, x) -1, (1.2)
where Q is a linear collision operator, and ρ(t, x) = R f (t, x, v)dv. The parameter α takes the value α = 1 in the diffusion scaling case and the value α = 0 in the high-field limit case. If periodic boundary conditions are considered, the function f satisfies the following boundary conditions (periodic in the spatial direction)

f (t, 0, v) = f (t, L, v), ∀v ∈ R, t ≥ 0.
In case of Dirichlet boundary conditions, we set

f (t, 0, v) = f L (t, v), ∀v ≥ 0, f (t, L, v) = f R (t, v), ∀v < 0, t ≥ 0, with lim v→±∞ f (t, x, v) = 0, ∀x ∈ [0, L], t ≥ 0.
In order to ensure the well-posedness of the problem, a boundary condition E(t, x = 0) has to be given or a mean electrostatic condition has to be added,

L 0 E(t, x)dx = λ(t), ∀t ≥ 0, (1.3)
where λ(t) is the potential drop at time t. Finally, an initial condition has to be considered for the transport equation

f (0, x, v) = f 0 (x, v), ∀x ∈ [0, L], v ∈ R. (1.4)
As said before, we are interested in two types of collision operators: the BGK operator for the diffusion limit (α = 1 in (1.1)):

Q(f )(t, x, v) = ρ(t, x) √ 2π exp(-v 2 /2) -f (t, x, v), with ρ(t, x) = R f (t, x, v)dv, (1.5)
and the Fokker-Planck operator for the high-field limit (α = 0 in (1.1))

Q(f )(t, x, v) = ∂ v [vf (t, x, v) + ∂ v f (t, x, v)]. (1.6)
We now recall that the Poisson equation (1.2) subjected to condition (1.3) implies the Ampère equation

∂ t E = - 1 ε α (J -J), (1.7)
where

J(t, x) = R vf (t, x, v)dv and J(t, x) = 1/L L 0 R vf (t, x, v)dvdx + λ ′ (t)/L, with λ defined in (1.3). Finally, integrating (1.2) between x = 0 and x = L gives E(t, 0) = E(t, L) ⇔ 1 L L 0 R f (t, x, v)dvdx = 1, ∀t ≥ 0, (1.8)
which translates the global neutrality of the plasma.

The main goal of this work is to extend the methodology of [START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF] and construct numerical schemes which are uniformly stable in the diffusion limit ε → 0 for system (1.1)-(1.2)-(1.5) with α = 1, and in the high-field limit for system (1.1)-(1.2)-(1.6) with α = 0. We will see that theses numerical schemes enjoy the following properties: (i) they are Asymptotic Preserving for the two asymptotics (Vlasov-Poisson-BGK towards (2.5) and Vlasov-Poisson-Fokker-Planck towards (3.11)), (ii) the high-field limit takes into account the first order correction terms (which means that the numerical scheme is first order (in ε) Asymptotic Preserving, (iii) they are not constrained by the space-diffusive CFL condition ∆t = O(∆x 2 ) for all regimes (ε << 1 and ε ≈ 1).

The rest of the paper is organized as follows: we first deal with the diffusion limit by presenting the derivation of the micro-macro model and its associated asymptotic preserving discretization. Next, we consider the high-field limit using the same outlines. Finally, the last part is devoted to the numerical results.

Diffusion limit

This section is devoted to the derivation and the discretization of the micro-macro model starting from the Vlasov-Poisson-BGK equation (1.1)-(1.2)-(1.5) with α = 1, which is well suited to study the diffusion limit. Obviously, the Fokker-Planck collision operator (1.6) can be used, but the BGK operator leads to simple computations.

We can define the null space of the linear BGK operator (1.5): using the notation f = R f (v)dv, it is given by N = Span{M } = {f = ρM, where ρ := f }, where M (v) = 1/ √ 2π exp(-v 2 /2) is the absolute Maxwellian and its rank R = (N ) ⊥ = {f such that f = 0}. Following [START_REF] Liu | Boltzmann equation: micromacro decompositions and positivity of shock profiles[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF], we decompose f as follows

f = ρM + g, with ρ(t, x) = R f (t, x, v)dv and M (v) = 1 √ 2π exp(-v 2 /2). (2.1)
As remarked in [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF], g = f -ρM is not necessarily small. We now introduce the orthogonal projector Π in L 2 (M -1 dv) onto N :

Πϕ = ϕ M, with ϕ = R ϕ(v)dv.
In the following, using these notations, the micro-macro decomposition model is derived and approximated to get an asymptotic preserving scheme in the diffusion limit.

Derivation of the micro-macro model

Introducing the transport operator T f = v∂ x f + E∂ v f , the decomposition (2.1) of f in the Vlasov equation (1.1) with α = 1 (diffusion scaling) leads to

∂ t (ρM ) + ∂ t g + 1 ε T (ρM ) + 1 ε T g = - 1 ε 2 g. (2.2)
Applying (I -Π) to (2.2) gives

(I -Π) ∂ t (ρM ) + 1 ε T (ρM ) + (I -Π) ∂ t g + 1 ε T g = - 1 ε 2 (I -Π)g.
Using the relations (I -Π)(∂ t (ρM )) = Π(g) = Π(∂ t g) = 0 (Lemma 3.1 of [START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF]), we get the microscopic part of the decomposition

∂ t g + 1 ε (I -Π)(T (ρM )) + 1 ε (I -Π)(T g) = - 1 ε 2 g.
Applying now Π to (2.2) leads to the macroscopic part of the decomposition

∂ t ρ + 1 ε ∂ x vg = 0. With the notations T f = v∂ x f + E∂ v f , M (v) = 1/ √ 2π exp(-v 2 /2
), the micro-macro formulation of (1.1)-(1.2) with α = 1 and Q chosen as the linear BGK operator (1.5) finally reads (1.5) with the initial data (1.4), then (ρ, g, E) = ( f , f -ρM, E) is a solution to (2.3) with the associated initial data

           ∂ t g + 1 ε (T g -∂ x vg M ) = 1 ε 2 [-g -εT (ρM )] , ∂ t ρ + 1 ε ∂ x vg = 0, ∂ x E = ρ -1.
2.1. (formal) (i) If (f, E) is a solution to (1.1)-(1.2)-
ρ 0 = f 0 , g 0 = f 0 -ρ 0 M.
(2.4)

(ii) Conversely, if (ρ, g, E) is a solution of (2.
3) with initial data ρ(t = 0) = ρ 0 and g(t = 0) = g 0 with g 0 = 0 then g(t) = 0, ∀t > 0 and (f = ρM + g, E) is a solution to (1.1)-(1.2)- (1.5) with initial data f 0 = ρ 0 M + g 0 .

Chapman-Enskog expansion

In this subsection, the formal derivation of the diffusion model is performed starting the micromacro model (2.3). Th diffusion model has been derived mathematically in [START_REF] Poupaud | Diffusion approximation of the linear semiconductor Boltzmann equation[END_REF]. Hereafter, we will see that the formal derivation is really straightforward starting from (2.3) (compared to the equivalent formulation of (1.1)-(1.2)-(1.5)), since the micro-macro model is well suited to deal with the asymptotic model in the diffusion limit. Indeed, for small ε, the first equation of (2.3) gives formally g = -εT

(ρM ) + O(ε 2 ) = -ε(vM ∂ x ρ + Eρ∂ v M ) + O(ε 2 ),
which, injected in the second equation of (2.3), leads to

∂ t ρ -∂ x v 2 M ∂ x ρ + Eρ v∂ v M = O(ε).
Since v 2 M = 1 and v∂ v M = -1, we finally get the asymptotic model (coupled with the Poisson equation)

∂ t ρ -∂ x (∂ x ρ -Eρ) = 0, ∂ x E = ρ -1.
(2.5)

Numerical scheme

We present here two discretizations of (2.3) which provide a numerical scheme to solve the original model (1.1)-(1.2)-(1.5) with α = 1 and have the following properties: (i) For all fixed ε > 0, the numerical schemes are consistent with (2.3) and (1.1)-(1.2)-(1.5); (ii) For fixed numerical parameters ∆t, ∆x, ∆v these schemes will degenerate into consistent discretizations of the asymptotic model (2.5) when ε → 0; moreover, for the second discretization, we have: (iii) The scheme is free from the usual diffusive CFL condition ∆t = O(∆x 2 ), and a time-implicit scheme is obtained for the diffusion term in (2.5) in the limit ε → 0.

Time explicit discretization

This subsection is devoted to the discretization of the micro-macro system (2.3) following the strategy proposed in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. We first present the time discretization and then tackle the full discretized problem.

Semi-discretization in time

The fixed time step is ∆t, and the current time t n = n∆t, n ∈ IN. Let ρ n and g n be the approximations of ρ and g at time t n . Using the notation

T f = v∂ x f +E n ∂ v f
, where E n satisfies the Poisson equation at time t n (∂ x E n = ρ n -1), the discretization of the first equation of (2.3) writes

g n+1 -g n ∆t + 1 ε (T g n -∂ x vg n M ) = 1 ε 2 -g n+1 -εT (ρ n M ) , (2.6)
with M (v) = 1/ √ 2π exp(-v 2 /2). Once g n+1 is computed thanks to (2.6), we look for the time approximation of the macroscopic part of (2.3). As in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF], we propose

ρ n+1 -ρ n ∆t + 1 ε ∂ x vg n+1 = 0. (2.7)
We now claim that, as ε → 0, the discretization (2.6)-(2.7) leads to a consistent time discretization of the asymptotic diffusion model (2.5). Indeed, when ε is small, we get from (2.6)

g n+1 = -εT (ρ n M ) + O(ε 2 ) = -ε(vM ∂ x ρ n + E n ρ n ∂ v M ) + O(ε 2 ).
Injecting it in (2.7) leads to

ρ n+1 -ρ n ∆t -∂ x (∂ x ρ n -E n ρ n ) = 0. (2.8)
Considering the Poisson equation ∂ x E n = ρ n -1, this corresponds to an explicit time discretization of the diffusion model (2.5).

Full-discretization

As in [START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF], we consider the micro equation of (2.3) at

x i+1/2 = (i+1/2)∆x, i = -1, • • • , N x and v k = v min + k∆v, k = 0, • • • , N v -1
where ∆x = L/N x and ∆v = (v maxv min )/N v (with v max = -v min ) are the phase space uniform discretization. Hence we denote by g n i+1/2,k an approximation of g(t n , x i+1/2 , v k ) and M k an approximation of M (v k ). On the other side, the macro part of (2.3) is approximated at

x i = i∆x, i = 0, • • • , N x so that we use the notations ρ n i ≈ ρ(t n , x i ) and ρ n i+1/2 = (ρ n i + ρ n i+1 )/2. The electric field is evaluated at x i+1/2 : E n i+1/2 ≈ E(t n , x i+1/2
) and the phase space fluxes are first order upwind fluxes

Φ i+1/2,k = v + k g n i+1/2,k -g n i-1/2,k ∆x + v - k g n i+3/2,k -g n i+1/2,k ∆x (2.9) Ψ i+1/2,k = E n,+ i+1/2 g n i+1/2,k -g n i+1/2,k-1 ∆v + E n,- i+1/2 g n i+1/2,k+1 -g n i+1/2,k ∆v , (2.10)
where the notations

v ± = (v ± |v|)/2 and E n,± i+1/2 = (E n i+1/2 ± |E n i+1/2 |)/2 are used.
Let us recall the given potential drop λ(t) given by (1.3) and the notation M k = ∆v Nv-1 k=0 M k . The following proposition presents the full discretized numerical scheme for the explicit version of the Asymptotic Preserving scheme in the diffusion limit.

Proposition 2.2. With the notations introduced above, the following discretization of the micromacro model (2.3) enjoys the Asymptotic Preserving property in the diffusion limit

g n+1 i+1/2,k ε 2 ∆t + 1 = g n i+1/2,k ε 2 ∆t -ε Φ i+1/2,k + Ψ i+1/2,k -Φ i+1/2,k M k -ε v k M k ρ n i+1 -ρ n i ∆x + E n i+1/2 ρ n i+1/2 M k+1 -M k-1 2∆v , (2.11)
ρ n+1 i -ρ n i ∆t = - 1 ε∆x v k (g n+1 i+1/2,k -g n+1 i-1/2,k ) , (2.12
)

E n i+1/2 -E n i-1/2 = (ρ n i -1) ∆x, with ∆x Nx i=0 E n i+1/2 = λ(t n ). ( 2

.13)

The associated asymptotic numerical scheme (ε → 0) is given by

ρ n+1 i -ρ n i ∆t - ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 + E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ∆x = 0. (2.14)
Together with (2.13), it provides a consistent discretization of the limit model (2.5). A necessary condition for the stability of (2.14) is ∆t = O(∆x 2 ).

Remark 2.1. For the micro part, the space and velocity derivatives for g n in (2.11) are approximated by a first order upwind scheme but many other choices of numerical fluxes are possible. The space derivative for ρ n and the velocity derivative for M in (2.11) are considered as a source term so that centered difference schemes are used. For the macro part, since it is approximated at point x i , the non-equilibrium flux can be naturally approximated by a central difference scheme. Obviously, other strategies are possible in which the same grid is used for the micro and the macro parts.

Proof. (Formal asymptotics) Let us now detail the formal proof of the asymptotic preserving property of the numerical scheme (2.11)-(2.12): we prove that the previous scheme degenerates into a consistent scheme for the diffusion equation (2.5) when ε → 0. In particular, we shall derive the asymptotic discretization of the diffusion model (2.5). From (2.11), we have for small ε

g n+1 i+1/2,k = -ε v k M k ρ n i+1 -ρ n i ∆x + E n i+1/2 ρ n i+1/2 M k+1 -M k-1 2∆v + O(ε 2 ),
which, injected in (2.12), leads to (up to first order in ε)

ρ n+1 i -ρ n i ∆t = - 1 ε∆x v k (g n+1 i+1/2,k -g n+1 i-1/2,k ) = 1 ∆x v 2 k M k ρ n i+1 -ρ n i ∆x - ρ n i+1 -ρ n i ∆x + 1 ∆x v k M k+1 -M k-1 2∆v (E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ) = ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 - E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2
∆x .

We used the following discrete integration by parts (neglecting the boundary terms)

v k M k+1 -M k-1 2∆v = v k-1 -v k+1 2∆v M k = -M k = -1.
The so-obtained discretization

ρ n+1 i -ρ n i ∆t - ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 + E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ∆x = 0, with ρ n i+1/2 = (ρ n i + ρ n i+1
)/2 and the discretization (2.13) of the Poisson equation is a consistent discretization of the limit model (2.5) where the diffusion term ∂ 2

x ρ is explicit.

Time implicit discretization

This part is devoted to the construction of a numerical scheme for the micro-macro model (2.3), the limiting scheme of which involves a time implicit treatment of the diffusion term ∂ 2 x ρ in (2.5). The motivation is to overcome the diffusion restriction on the time step ∆t = O(∆x 2 ).

Proposition 2.3. With the notations introduced in Proposition 2.2 and (2.13), the following discretization of the micro-macro model (2.3) enjoys the Asymptotic Preserving property and it is free from the usual diffusion condition on the time step

g n+1 i+1/2,k ε 2 ∆t + 1 = g n i+1/2,k ε 2 ∆t -ε Φ i+1/2,k + Ψ i+1/2,k -Φ i+1/2,k M k -ε v k M k ρ n i+1 -ρ n i ∆x + E n i+1/2 ρ n i+1/2 M k+1 -M k-1 2∆v , (2.15)
ρ n+1 i -ρ n i ∆t = - 1 ε∆x 1 + ∆t ε 2 -1 v k (g n i+1/2,k -g n i-1/2,k ) - ∆t ε v k (Φ i+1/2,k -Φ i-1/2,k ) - ∆t ε v k (Ψ i+1/2,k -Ψ i-1/2,k ) - ∆t ε v 2 k M k ρ n+1 i+1 -2ρ n+1 i + ρ n+1 i-1 ∆x - ∆t ε (E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ) . (2.

16)

The associated asymptotic numerical scheme (ε → 0) is given by

ρ n+1 i -ρ n i ∆t - ρ n+1 i+1 -2ρ n+1 i + ρ n+1 i-1 ∆x 2 + E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ∆x = 0, (2.17)
which, together with (2.13), provides a consistent discretization of the limit model (2.5).

Proof. To derive this implicit discretization of (2.3), let us start with the semi-discrete case. First of all, the origin of the diffusion term has to be identified. Following the Chapman-Enskog procedure detailed in subsection 2.2, for small ε the kinetic part g is given by

g = -εT (ρM ) = -ε(vM ∂ x ρ + Eρ∂ v M ) + O(ε 2 ).
Therefore (-εvM ∂ x ρ) is the term that generates the diffusion ∂ 2 x ρ in the asymptotics ε → 0. In the following, we exploit this to derive an implicit treatment of the spatial diffusion.

First, we rewrite (2.6) to obtain the expression of g n+1 (recalling that T g n = v∂ x g n + E n ∂ v g n ),

g n+1 = 1 + ∆t ε 2 -1 g n - ∆t ε (T g n -∂ x vg n M ) - ∆t ε vM ∂ x ρ n - ∆t ε E n ρ n ∂ v M . (2.18)
We can inject this expression in the flux of the macro equation to get

ρ n+1 -ρ n ∆t + 1 ε 1 + ∆t ε 2 -1 ∂ x v g n - ∆t ε (T g n -∂ x vg n M ) - ∆t ε vM ∂ x ρ n - ∆t ε E n ρ n ∂ v M = 0.
(2.19) Some computations enable to simplify the term into brackets

v g n - ∆t ε (T g n -∂ x vg n M ) - ∆t ε vM ∂ x ρ n - ∆t ε E n ρ n ∂ v M = vg n - ∆t ε v 2 ∂ x g n + vE n ∂ v g n + ∆t ε vM ∂ x vg n - ∆t ε v 2 M ∂ x ρ n - ∆t ε v∂ v M E n ρ n = vg n - ∆t ε ∂ x v 2 g n - ∆t ε ∂ x ρ n + ∆t ε E n ρ n ,
where some terms vanish (third and fourth terms of the second line) using zero-mean average property for g n and symmetry property for M in the velocity direction. Then the numerical scheme for the macro part ρ can be written as

ρ n+1 -ρ n ∆t + 1 ε 1 + ∆t ε 2 -1 ∂ x vg n - ∆t ε ∂ x v 2 g n - ∆t ε ∂ x ρ n + ∆t ε E n ρ n .
Then, it is possible to make implicit the diffusion term to finally get

ρ n+1 -ρ n ∆t + 1 ε 1 + ∆t ε 2 -1 ∂ x vg n - ∆t ε ∂ x v 2 g n - ∆t ε ∂ x ρ n+1 + ∆t ε E n ρ n . (2.20)
The full discretization of this new scheme can be done as previously. Indeed, starting with (2.11) which we inject into (2.12), the full discretization of the macro part becomes

ρ n+1 i -ρ n i ∆t = - 1 ε∆x 1 + ∆t ε 2 -1 v k (g n i+1/2,k -g n i-1/2,k ) - ∆t ε v k (Φ i+1/2,k -Φ i-1/2,k ) - ∆t ε v k (Ψ i+1/2,k -Ψ i-1/2,k ) - ∆t ε v 2 k M k ρ n i+1 -2ρ n i + ρ n i-1 ∆x - ∆t ε (E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ) ,
where the fluxes Φ i+1/2,k and Ψ i+1/2,k are given by (2.9) and (2.10). Then, it is possible to make implicit the diffusion term to finally get

ρ n+1 i -ρ n i ∆t = - 1 ε∆x 1 + ∆t ε 2 -1 v k (g n i+1/2,k -g n i-1/2,k ) - ∆t ε v k (Φ i+1/2,k -Φ i-1/2,k ) - ∆t ε v k (Ψ i+1/2,k -Ψ i-1/2,k ) - ∆t ε v 2 k M k ρ n+1 i+1 -2ρ n+1 i + ρ n+1 i-1 ∆x - ∆t ε (E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ) , (2.21)
The study of the asymptotic ε → 0 is similar to the explicit case. As ε → 0, g n i+1/2,k = O(ε) so that the fluxes Ψ i+1/2,k and Φ i+1/2,k are also of order ε. Moreover, (1 + ∆t/ε 2 ) -1 is equivalent to ε 2 /∆t when ε << 1, so that, at the limit ε = 0, (2.21) reduces to

ρ n+1 i -ρ n i ∆t - ρ n+1 i+1 -2ρ n+1 i + ρ n+1 i-1 ∆x 2 + E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ∆x = 0,
which, together with (2.13) is a consistent discretization of the limit model (2.5); the diffusion term has been implicited so that this scheme (2.15)-(2.16) is free from the standard diffusion condition on the time step.

Boundary conditions

In this subsection, boundary conditions in the x variable are investigated. As explained in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF], the original boundary conditions on f cannot always be translated into boundary conditions on the macro part ρ and the micro part g separately. For example, while periodic boundary conditions are straightforwardly transposed into boundary conditions on ρ and g, Dirichlet boundary conditions cannot be extracted for ρ and g and a special handling is needed in this case. This is because only the ingoing part of f can be imposed as a Dirichlet boundary condition, as given by (1.3).

For the full kinetic model (1.1), the boundary conditions we consider here are the following (for given functions f L and f R )

f (t, x = 0, v) = f L (v), v > 0, f (t, x = L, v) = f R (v), v < 0,
which has the following discrete version (x 0 = 0 and

x Nx = L) f (t, x 0 , v k ) = f L (v k ), v k > 0, f (t, x Nx , v k ), v k < 0.
For the micro-macro decomposition, it becomes for all t ≥ 0

ρ(t, x 0 )M k + 1 2 (g(t, x 1/2 , v k ) + g(t, x -1/2 , v k )) = f L (v k ), v k > 0, (2.22) ρ(t, x Nx )M k + 1 2 (g(t, x Nx+1/2 , v k ) + g(t, x Nx-1/2 , v k )) = f R (v k ), v k < 0. (2.23)
For the other velocities, we impose Neumann conditions which can be viewed at a first order approximation of g at the boundary

g(x -1/2 , v k ) = g(x 1/2 , v k ), v k < 0, (2.24) g(x Nx+1/2 , v k ) = g(x Nx-1/2 , v k ), v k > 0. (2.25)
Hence, following the procedure introduced in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF], we write the discrete macro equation (2.12) for i = 0 which involves g n+1 -1/2,k for all k, and use (2.22) and (2.24) at t = t n+1 to get

ρ n+1 0 -ρ n 0 ∆t + 1 ∆x v k g n+1 1/2,k -v + k [2(f L (v k ) -ρ n+1 0 M k ) -g n+1 1/2,k ] -v - k g n+1 1/2,k = 0. (2.26)
In the same spirit for i = N x , we get using (2.23) and (2.25

) at t = t n+1 ρ n+1 Nx -ρ n Nx ∆t + 1 ∆x v + k g n+1 Nx-1/2,k + v - k [2(f R (v k ) -ρ n+1 Nx M k ) -g n+1 Nx-1/2,k ] -v k g n+1 Nx-1/2,k = 0. (2.27)
Finally, using the boundary conditions (2.22) and (2.24), we can update

g n+1 -1/2,k for all k g n+1 -1/2,k = g n+1 1/2,k , v k < 0, g n+1 -1/2,k = 2(f L (v k ) -ρ n+1 0 M k ) -g n+1 1/2,k , v k > 0, (2.28)
and g n+1 Nx+1/2,k using (2.23) and (2.25)

g n+1 Nx+1/2,k = g n+1 Nx-1/2,k , v k > 0, g n+1 Nx+1/2,k = 2(f R (v k ) -ρ n+1 Nx M k ) -g n+1 Nx-1/2,k , v k < 0. (2.29)
Finally, assuming g n i+1/2,k is known for i = -1, • • • , N x and ρ n i is known for i = 0, • • • , N x , the procedure can be summarized into 4 points

• advance (2.11) for i = 0, • • • , N x -1 • advance (2.12) for i = 1, • • • , N x -1 • compute ρ n+1 0 and ρ n+1
Nx with (2.26) and (2.27)

• compute g n+1 -1/2,k and g n+1 Nx+1/2,k with (2.28) and (2.29). For the Poisson equation, the boundary conditions of Dirichlet type are imposed on the electric potential so that (2.13) is used for the computation of the electric field.

Remark 2.2. Note that this approach can also be applied to solve the micro-macro system where the macro part (2.16) is now discretized with an implicit scheme.

The procedure is similar to the one presented just above. Once the micro part is advanced using (2.11)- (2.16) for i = 0, • • • , N x -1, we can compute ρ n+1 at the boundaries using (2.26) and (2.27). Then, the system coming from (2.16) to be solved is the following, for i = 2, . . . , N x -2

(Aρ n+1 ) i = ρ n i - ∆t ε∆x v k (g n+1 i+1/2,k -gn+1 i-1/2,k ) , i = 2, . . . , N x -2 (2.

30)

where A is the tridiagonal matrix A = Tridiag(-α, 1 + 2α, -α) with α = -∆t 2 /(∆x 2 (ε 2 + ∆t 2 )) and

gn+1 i+1/2,k = g n+1 i+1/2,k + ε ε 2 /∆t + 1 v k M k ρ n i+1 -ρ n i ∆x ,
with g n+1 i+1/2,k is computed with (2.11). The system is tridiagonal since ρ n+1 0 and ρ n+1 Nx have been computed using (2.26) and (2.27); hence the first line of the right hand side of the linear system (2.30) is given by

(1 + 2α)ρ n+1 1 -αρ n+1 2 = ρ n 1 - ∆t ε∆x v k (g n+1 3/2,k -gn+1 1/2,k ) + αρ n+1 0 ,
whereas the last line is given by

-αρ n+1 Nx-2 + (1 + 2α)ρ n+1 Nx-1 = ρ n Nx-1 - ∆t ε∆x v k (g n+1 Nx-1/2,k -gn+1 Nx-3/2,k ) + αρ n+1 Nx .
Remark 2.3. Note that the staggered grids are not necessary. They provided a simplified framework but the same grid can be adopted for the micro and the macro part of the model ; this may be interesting for multidimensional applications.

High-field limit

This section is devoted to the derivation of the micro-macro model in the high-field asymptotic limit, and of asymptotic preserving schemes in this context. The starting model is the Vlasov-Poisson-Fokker-Planck equation (1.1)-(1.2)-(1.6) with α = 0 which we re-write as follows

∂ t f + v∂ x f = 1 ε ∂ v [(v -E)f + ∂ v f ] =: 1 ε Lf. (3.1)
Considering the linear operator L in the right hand side of (3.1), we can define its null space which is given by N (L) = Span{M} = {f = ρM, where ρ := f }, where

M(v) = 1/ √ 2π exp(-(v - E) 2 /2) is the shifted Maxwellian and its rank R(L) = (N ) ⊥ (L) = {f such that f = R f dv = 0}.
Following [START_REF] Liu | Boltzmann equation: micromacro decompositions and positivity of shock profiles[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF], we decompose f as follows

f = ρM + g, (3.2)
with ρ(t, x) = R f (t, x, v)dv and M(t, x, v) = 1 √ 2π exp(-(v -E(t, x)) 2 /2). (3.3)
As in the diffusion case, we introduce the orthogonal projector Π in L 2 (M -1 dv) onto N (L):

Πϕ = ϕ M, with ϕ = R ϕ(v)dv. (3.4)
Note that here M is not an absolute Maxwellian as previously, but depends on (t, x) throught the electric field E. The rest of this section is devoted to the derivation and discretization of the micro-macro decomposition model corresponding to the high-field scaling.

Derivation of the micro-macro model

Introducing the micro-macro decomposition (3.2) in the Vlasov equation (3.1), we obtain

∂ t (ρM) + ∂ t g + v∂ x (ρM) + v∂ x g = 1 ε Lg. (3.5)
Applying (I -Π) to (3.5) yields

(I -Π)[∂ t (ρM) + v∂ x (ρM)] + (I -Π)[∂ t g + v∂ x g] = 1 ε (I -Π)Lg. Since Π(g) = Π(∂ t g) = Π(Lg) = 0 (Lemma 3.1 of [3]
), and

Π[∂ t (ρM)] = Π [M∂ t ρ + ρ(∂ t E)(v -E)M] = M∂ t ρ,
we get, using the Ampère equation (1.7) (we recall the notation J = 1/L L 0 J(t, x)dx + λ ′ (t)/L with λ defined in (1.3))

∂ t g + (I -Π)(v∂ x g) = 1 ε Lg -ε(I -Π)[v∂ x (ρM)] + ε(J -J)(v -E)(ρM) . (3.6)
The current J can be defined as a function of ρ, E and g using the following expression

J(t, x) = R v [ρ(t, x)M(t, x, v) + g(t, x, v)] dv = ρ(t, x)E(t, x) + R vg(t, x, v)dv.
Applying now the Π projection to (3.5) and using its expression (3.4), we can write the micro-macro formulation of (1.1)-(1.2)-(1.6) in the high-field case (1.6) with the initial data (1.4), then (ρ, g, E) = ( f , f -ρM, E) is a solution of (3.7) with the associated initial data

       ∂ t g + (v∂ x g -∂ x vg M) = 1 ε Lg -εv∂ x (ρM) + ε∂ x v(ρM) M + ε(J -J)(v -E)(ρM) , ∂ t ρ + ∂ x vM + ∂ x vg = 0, ∂ x E = ρ -1. ( 3
) (i) If (f, E) is a solution of (1.1)-(1.2)-
ρ 0 = f 0 , g 0 = f 0 -ρ 0 M(t = 0). (3.8) (ii) Conversely, if (ρ, g, E)
is a solution of (3.7) with initial data ρ(t = 0) = ρ 0 and g(t = 0) = g 0 with g 0 = 0 then g(t) = 0, ∀t > 0 and

(f = ρM + g, E) is a solution to (1.1)-(1.2)-(1.6) with initial data f 0 = ρ 0 M(t = 0) + g 0 .

Chapman-Enskog expansion

In the same spirit as previously, we shall derive formally the asymptotic high-field model. Let us refer to [START_REF] Bostan | Electric high-electric field limit for the Vlasov-Maxwell-Fokker-Planck system[END_REF][START_REF] Nieto | High-field limit for the Vlasov-Poisson-Fokker-Planck system[END_REF] for a mathematical derivation. Starting from (3.7) which is an equivalent formulation of (1.1)-(1.2)-(1.6), we will see that the asymptotic high-field model can be quite easily obtained since the micro-macro model (3.7) is a well-suited model to do that. Indeed, for small ε, the first equation of (3.7) gives formally

g = L -1 εv∂ x (ρM) -εM∂ x v(ρM) -ε(J -J)(v -E)(ρM) + O(ε 2 ) = εL -1 [vM∂ x ρ + ρv∂ x E(v -E)M -M∂ x (ρE) -(J -J)(v -E)ρM] + O(ε 2 ) = ε(∂ x ρ -ρ(J -J)) L -1 ((v -E)M) -ερ∂ x E L -1 ((1 -v(v -E))M) + O(ε 2 ). (3.9)
where L -1 is the pseudo-inverse of L (i.e. the inverse when restricted to N ⊥ (L) = R(L)). This has to be injected in the second equation of (3.7).

To get an explicit expression of the pseudo-inverse L -1 , let us compute L(a(v)M) for a given function a depending only on v

L(a(v)M) = ∂ v [(v -E)aM + ∂ v (aM)] = ∂ v [(v -E)aM + a ′ (v)M -a(v)(v -E)M] = ∂ v [a ′ (v)M] = a ′′ (v)M -a ′ (v)(v -E)M.
For a(v) = v, v 2 /2, we then obtain the following relations (using

Lg = f =⇒ g -g M = L -1 f since L is invertible on N ⊥ (L) = R(L) = {f such that f = R f dv = 0}) L(vM) = -(v -E)M then -vM + vM M = L -1 ((v -E)M) and L v 2 2 M = M -v(v -E)M then v 2 2 M - v 2 2 M M = L -1 ((1 -v(v -E))M).
Then, injecting these two last expressions into (3.9) leads to

g = -ε(∂ x ρ -ρ(J -J))(v -E)M -ερ∂ x E v 2 2 M -M (1 + E 2 ) 2 + O(ε 2 ) = -ε(∂ x ρ -ρ(J -J))(v -E)M - ε 2 ρ∂ x E (v 2 -1 -E 2 )M + O(ε 2 ). ( 3.10)
We are now ready to write the asymptotic model using the previous computations. We first use (3.10) in the second equation of (3.7) and get

∂ t ρ + ∂ x v(ρM) + ε∂ x v -(∂ x ρ -ρ(J -J))(v -E)M -ρ∂ x E(v 2 /2 -(1 + E 2 )/2)M = O(ε 2 ).
The third term of this equation gives

-ε∂ x (∂ x ρ -ρ(J -J )) v(v -E)M = -ε∂ 2 x ρ + ε∂ x (ρ(J -J)),
whereas the fourth one writes

-ε∂ x ρ∂ x E v 3 2 M - 1 + E 2 2 vM = -ε∂ x ρ∂ x E 1 2 (3E + E 3 -E -E 3 ) = -ε∂ x (ρE∂ x E).
Near equilibrium, the current J can be expressed as a function of ρ and E only ; indeed, since from (3.10), g = O(ε) and we have

J = R vf dv = ρ R vMdv + O(ε) = ρE + O(ε).
With the conditions (1.8) and (1.3), and using the Poisson equation, we have :

J = 1 L L 0 ρEdx + 1 L λ ′ (t) + O(ε) = 1 L L 0 (∂ x E + 1)Edx + 1 L λ ′ (t) + O(ε) = 1 2L L 0 ∂ x (E 2 )dx + 1 L L 0 Edx + 1 L λ ′ (t) + O(ε) = 1 L (λ(t) + λ ′ (t)) + O(ε).
Finally, using the Poisson equation ∂ x E = ρ -1, the asymptotic model writes (with λ given by

(1.3))    ∂ t ρ + ∂ x (ρE) = ε∂ 2 x ρ -ε∂ x (ρE) + ε λ + λ ′ L ∂ x ρ, ∂ x E = ρ -1. (3.11)
If the condition λ(t) = 0 is fulfilled, we recover the standard high-field model derived in [START_REF] Bonilla | High-field limit of the Vlasov-Poisson-Fokker-Planck system for different perturbation methods[END_REF].

Numerical scheme

We then present a discretization of (3.7) which provides a numerical scheme to solve the original model (1.1)-(1.2)-(1.6) with α = 0 and has the following properties: (i) For all fixed ε > 0, the numerical scheme is consistent with (3.7) and (1.1)-(1.2)-(1.6); (ii) For fixed numerical parameters ∆t, ∆x, ∆v the scheme will degenerate into a consistent discretization of the asymptotic model (3.11) when ε → 0, up to terms of order O(ε 2 ). We consider the same outlines as in the diffusion case (subsection 2.3).

Time explicit discretization

This subsection is devoted to the discretization of the micro-macro system (3.7) following the strategy proposed in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. We first present the time discretization and then tackle the full discretized problem.

Semi-discretization in time

As previously, we denote by ∆t the time step, t n = n∆t, n ∈ IN the current time, g n , ρ n and E n the approximations of g, ρ and E at time t n . The semi-implicit time discretization of the micro part of (3.7) is

g n+1 -g n ∆t + v∂ x g n -∂ x vg n M n = 1 ε Lg n+1 -εv∂ x (ρ n M n ) + ε∂ x v(ρ n M n ) M n + ε(J n -Jn )(v -E n )(ρ n M n ) (3.12) with M n = 1/ √ 2π exp(-(v -E n ) 2 /2).
The macro part of (3.7) is discretized as follows

ρ n+1 -ρ n ∆t + ∂ x v(ρ n M n ) + ∂ x vg n+1 = 0. (3.13)
We now prove formally that the semi-discretization (3.12)-(3.13) leads to a consistent time discretization of the limit model (3.11) when ε → 0, ∆t being fixed. Indeed, from (3.12), we get when ε is small,

g n+1 = L -1 εv∂ x (ρ n M n ) -ε∂ x v(ρ n M n ) M n -ε(J n -Jn )(v -E n )(ρ n M n ) + O(ε 2 ).
We then express (3.13) in terms of ρ only by injecting the previous equality in (3.13)

ρ n+1 -ρ n ∆t + ∂ x v(ρ n M n ) + ε∂ x vL -1 v∂ x (ρ n M n ) -∂ x v(ρ n M n ) M n -(J n -Jn )(v -E n )(ρ n M n ) = O(ε 2 ).
Using computations which are similar to the continuous case, we finally get the asymptotic time discretization model

ρ n+1 -ρ n ∆t + ∂ x (E n ρ n ) = ε∂ 2 x ρ n -ε∂ x (ρ n E n ) .
Considering the Poisson equation ∂ x E n = ρ n -1, this is an explicit time discretization of the continuous asymptotic model (3.11).

Full-discretization As in the diffusion case, we consider the micro equation of (3.7) at

x i+1/2 = (i + 1/2)∆x, i = -1, • • • , N x and v k = v min + k∆v, k = 0, • • • , N v -1
where ∆x = L/N x and ∆v = (v maxv min )/N v (with v max = -v min ) are the phase space uniform discretization. Hence we denote by g n i+1/2,k an approximation of g(t n , x i+1/2 , v k ). On the other side, the macro part (3.7) is approximated at

x i = x min + i∆x, i = 0, • • • , N x and ρ n
i denotes an approximation of ρ(t n , x i ). The electric field is evaluated at x i+1/2 : E n i+1/2 ≈ E(t n , x i+1/2 ) as well as the current J n i+1/2 ≈ J(t n , x i+1/2 ) whereas for its spatial average, we have Jn ≈ Jn ; the Maxwellian is approximated as follows

M n i+1/2,k = 1/ √ 2π exp(-(v k -E n i+1/2 ) 2 /2). It enables to determine L i+1/2
such that for a function on the velocity grid (g k ) k , we have

L i+1/2 g k ≈ (Lg)(x i+1/2 , v k ). Finally, for the macroscopic flux ∂ x v(ρ n M n ) = ∂ x (ρ n E n ),
we use a Deshpande-Pullin [START_REF] Mandal | Kinetic flux vector splitting for Euler equations[END_REF] (but any other more sophisticated choice of numerical fluxes can be simply adapted)

F i+1/2 (ρ n ) = v + k ρ n i M n i,k + v - k ρ n i+1 M n i+1,k . (3.14)
Proposition 3.2. With the notations introduced above, the following discretization of the micromacro model (3.12) enjoys the Asymptotic Preserving property in the high-field limit

1 - ∆t ε L i+1/2 g n+1 i+1/2,k = g n i+1/2,k -∆t Φ i+1/2,k -Φ i+1/2,k M n i+1/2,k + M n i+1/2,k (v k -E n i+1/2 ) ρ n i+1 -ρ n i ∆x + G n i+1/2,k (3.15) ρ n+1 -ρ n ∆t + F i+1/2 (ρ n ) -F i-1/2 (ρ n ) ∆x + 1 ∆x v k (g n+1 i+1/2,k -g n+1 i-1/2,k ) = 0, (3.16) E n i+1/2 -E n i-1/2 = (ρ n i -1) ∆x, with ∆x Nx i=0 E n i+1/2 = λ(t n ), (3.17)
where the space flux Φ i+1/2,k is given by (2.9), the macroscopic flux F i+1/2 (ρ n ) is given by (3.14) and

G i+1/2,k = ρ n i+1/2 M n i+1/2,k ((E n i+1 -E n i )/∆x [v k (v k -E n i+1/2 ) -1] -(J n i+1/2 -Jn )(v k -E n i+1/2
)). The associated asymptotic numerical scheme (ε → 0) is given by

ρ n+1 i -ρ n i ∆t + E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2 ∆x = ε ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 -ε E n i+1/2 ρ n i+1/2 -E n i-1/2 ρ n i-1/2
∆x , (3.18) Together with (3.17), it is a consistent discretization of the limit model (3.11). A necessary condition for the stability of (3.18)

is ∆t = O(∆x 2 /ε).
Proof. Before proving (formally) that the numerical scheme (3.15)-(3.16)-(3.17) is asymptotic preserving, let us detail the discretization of the collision operator L.

For the discretization of Lg, we apply the strategy of [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 in the Coulomb case[END_REF]: the Fokker-Planck operator (3.1) is rewritten as

(Lg)(v) = ∂ v [M∂ v (g/M)], with M(v) = 1 √ 2π exp(-(v -E) 2 /2).
Starting with this formulation, we employ the following velocity discretization

(Lg) k = 1 2 D + MD -(g/M) k k + 1 2 D -MD + (g/M) k k , (3.19)
where

D + (g) k = (g k+1 -g k )/∆v and D -(g) k = (g k -g k-1
)/∆v are discrete differential operators.

The expression (3.19) can also be transformed into

(Lg)(v k ) ≈ (Lg) k = 1 2∆v 2 [g k-1 (1 + M k /M k-1 ) + g k (-2 -(M k+1 + M k-1 )/M k ) + g k+1 (1 + M k /M k+1 )] . (3.20)
We then approximate (Lg) k by a matrix vector product (Ag) k where A is a N v × N v tridiagonal matrix and g the vector of components g k , k = 0, ..., N v -1. The pseudo-inverse L -1 f of g is obtained from the following operation

f = Lg ⇐⇒ L -1 f = g -g M,
so that, at the discrete level, we have

f k = (Ag) k ⇐⇒ (A -1 f ) k = g k - Nv-1 k=0 g k ∆v M k .
Now, let us prove (formally) that the numerical scheme (3.15)-(3.16)-(3.17) enjoys the asymptotic preserving property. To do that, we first observe that (1 -∆t/ε L) -1 = -ε/∆t L -1 + O(ε 2 ), so that (3.15) becomes, for small ε

g n+1 i+1/2,k = εL -1 i+1/2 v k M n i+1/2,k ρ n i+1 -ρ n i ∆x + ρ n i+1/2 v k E n i+1 -E n i ∆x (v k -E n i+1/2 )M n i+1/2,k -M n i+1/2,k E n i+1/2 ρ n i+1 -ρ n i ∆x -M n i+1/2,k ρ n i+1/2 E n i+1 -E n i ∆x -ρ n i+1/2 (J n i+1/2 -Jn )(v k -E n i+1/2 )M n i+1/2,k = εL -1 i+1/2 (v k -E n i+1/2 )M n i+1/2,k ρ n i+1 -ρ n i ∆x -ρ n i+1/2 (J n i+1/2 -Jn ) + εL -1 i+1/2 (v 2 k -v k E n i+1/2 -1)M n i+1/2,k ρ n i+1/2 E n i+1 -E n i ∆x .
Considering the inversion of L i+1/2 as in the continuous case (i.e. exact with respect to v), we get

g n+1 i+1/2,k = -ε(v k -E n i+1/2 )M n i+1/2,k ρ n i+1 -ρ n i ∆x -ρ n i+1/2 (J n i+1/2 -Jn ) - ε 2 (1 + (E n i+1/2 ) 2 -v 2 k )M n i+1/2,k ρ n i+1/2 E n i+1 -E n i ∆x .
We have now to inject this expression in the macro part, i.e. we have to compute v k (g n+1 i+1/2,k -

g n+1 i-1/2,k ) . First we compute v k g n+1 i+1/2,k = ∆v Nv k=0 (v k g n+1 i+1/2,k ) v k g n+1 i+1/2,k = -ε ρ n i+1 -ρ n i ∆x -ρ n i+1/2 (J n i+1/2 -Jn ) v k (v k -E n i+1/2 )M n i+1/2,k - ε 2 ρ n i+1/2 E n i+1 -E n i ∆x v k 1 + (E n i+1/2 ) 2 -v 2 k M n i+1/2,k = -ε ρ n i+1 -ρ n i ∆x -ρ n i+1/2 (J n i+1/2 -Jn ) -εE n i+1/2 ρ n i+1/2 E n i+1 -E n i ∆x .
For the last equality, we used the following identities

v k (v k -E n i+1/2 )M n i+1/2,k = (v k -E n i+1/2 ) 2 + E n i+1/2 (v k -E n i+1/2 ) M n i+1/2,k = 1,
and

v k 1 + (E n i+1/2 ) 2 -v 2 k M n i+1/2,k = v k M n i+1/2,k + (E n i+1/2 ) 2 v k M n i+1/2,k -v 3 k M n i+1/2,k = -2E n i+1/2 .
Now, as in the continuous case (when λ = 0 for simplicity), we remark that

J n i+1/2 -Jn = ρ n i+1/2 E n i+1/2 + O(ε), so we can replace (J n i+1/2 -Jn ) by ρ n i+1/2 E n i+1/2 since the correction terms are of order O(ε 2 ). Considering the term v k (g n+1 i+1/2,k -g n+1 i-1/2,k ) enables to conclude 1 ∆x v k (g n+1 i+1/2,k -g n+1 i-1/2,k ) = -ε ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 + ε (ρ n i+1/2 ) 2 E n i+1/2 -(ρ n i-1/2 ) 2 E n i-1/2 /∆x - ε ∆x 2 ρ n i+1/2 (E n i+1 -E n i )E n i+1/2 -ρ n i-1/2 (E n i -E n i-1 )E n i-1/2 = -ε ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 + ε ∆x   ρ n i+1/2 E n i+1/2 =1 using Poisson equation (ρ n i+1/2 -(E n i+1 -E n i )/∆x) -ρ n i-1/2 E n i-1/2 (ρ n i-1/2 -(E n i -E n i-1 )∆x) =1 using Poisson equation     = -ε ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 + ε ∆x ρ n i+1/2 E n i+1/2 -ρ n i-1/2 E n i-1/2 .
This is a discretization of first order terms in the limit model (3.11), which proves that our scheme is an asymptotic preserving scheme and recover also the first order corrective terms in ε.

Remark 3.1. Let us remark that a time implicit discretization is also possible in this case, as for the diffusion asymptotic. The approach is similar (but more technical) as for the diffusion case: the diffusion term ε∂ 2 x ρ is provided by vM∂ x ρ and by ME∂ x ρ. In fact, this procedure is not necessary for small ε since the diffusive CFL condition for the asymptotic model is of the form ∆t = O(∆x 2 /ε), which is not restrictive in this case. We observed that the implicit or explicit approach give rise to equivalent results. The implicit approach is however used in the numerical results (it ensures that the time step only is only restricted by the transport condition).

Numerical results

We will test our algorithms with the classical test cases of plasma physics.

Periodic boundary condition and Landau damping

The initial condition associated to the scaled Vlasov-Poisson-BGK equation has the following form

f 0 (x, v) = 1 √ 2π exp(-v 2 /2)(1 + α cos(kx)), (x, v) ∈ [0, 2π/k] × R,
where k = 0.5 is the wave number and α = 0.05 is the amplitude of the perturbation. A cartesian mesh is used to represent the phase space with a computational domain [0

, 2π/k] × [v min , v max ], -v min = v max = 6.
The number of mesh points in the spatial and velocity directions is designated by N x and N v respectively.

Diffusion case Our purpose is to compare the three models considering the diffusion scaling: (i) Vlasov refers to the Vlasov-Poisson-BGK equation (1.1)-(1.2)-(1.5) with α = 1, discretized in an explicit way (see [START_REF] Crouseilles | A hybrid kinetic-fluid model for solving the Vlasov-BGK equations[END_REF] for instance), (ii) MM stands for the micro-macro associated model (2.3) discretized using the implicit discretization of Proposition 2.3 and (iii) LIM refers to the diffusion model (2.5) discretized in an explicit way (2.14). For these three models, the time step is chosen to satisfy the CFL conditions (induced by the transport and the diffusion); using the notation β x = ∆x/v max and β v = ∆v/E max (where E max denotes the maximum value of the electric field) we choose

-Vlasov: ∆t = C min εβ x , εβ v , ε 2 , -MM: ∆t = min (Cεβ x / max(0, 1 -Cβ x /ε), Cεβ v / max(0, 1 -Cβ v /ε)), -LIM: ∆t = C min ∆x/E max , ∆x 2 .
The CFL number C is chosen equal to C = 0.5. Finally, the phase space grid is discretized using N x = N v = 128 points except for ε = 1 in which 512 points are used per direction.

Let us remark that the use of the time implicit discretization (given by Proposition 2.3) for the micro-macro model gives rise to numerical results which are very close to those produced by the explicit discretization (given by Proposition 2.2). Hence, only the implicit discretization is used for the following results.

For our micro-macro model, the initial conditions are given by ρ(t = 0, x) = 1 + α cos(kx), and g(t = 0, x, v) = 0, whereas for the diffusion model the initial condition is

ρ(t = 0, x) = 1 + α cos(kx).
We are looking for the time evolution of the electric energy E(t) L 2 in log scale for the three models (Vlasov, MM and LIM) but we also plot the density ρ and the electric field E as a function of x for different times. We are interested in the influence of the ε parameter. On the one hand, when ε ≈ 1, the MM model is very close to the reference Vlasov solution (see Figures 2 and3). The same is true for ε = 0.5 (see Figure 4). On the other hand, with a fixed ∆t, ∆x and ∆v, we observe that the MM model gives the right limit ε → 0 since its solution is very close to the LIM model. Indeed, on Figures 5 and6, we can observe that even for the density and the electric field, both models (MM and LIM) give very similar results. Obviously, when ε is small, the computational time required for the explicit Vlasov scheme is too large. For some given numerical parameters, the CPU time of the micro-macro model is about twice the CPU time of the simulation of the Vlasov-Poisson-BGK equation. If ε is very small, the explicit computation (Vlasov) becomes unreachable while the MM scheme coincides with the LIM model with nongrowing computational cost. With the same geometry as previously, the numerical parameters are chosen as follows:

N x = N v = 256.
For our micro-macro model (3.7), the initial conditions are given by ρ(t = 0, x) = 1 + α cos(kx), and g(t = 0, x, v) = f (t = 0, x, v)ρ(t = 0, x)M(t = 0, x, v),

with M(t = 0, x, v) = 1/ √ 2π exp(-(v -E(t = 0, x)) 2 /2
) where E(t = 0, x) is the electric field satisfying the Poisson equation at initial time ∂ x E(t = 0, x) = ρ(t = 0, x) -1.

We are looking for the time evolution of the electric energy E(t) L 2 in log scale for the three models (Vlasov, MM and LIM). We are interested in the influence of the ε parameter. On the one side, when ε = 1, 0.5, the MM model is very close to the reference Vlasov solution when regarding the electric energy, the density or the electric field (see Figures 7 and8 for ε = 1 and Figures 9 and 10 for ε = 0.5). On the other side, with a fixed ∆t, ∆x and ∆v, we observe that the MM model gives the right limit ε → 0 (see Figures 11 and12 for ε = 0.01), since its solution is very close to the LIM model. As claimed before, the computational time is independent from ε since the time step and the mesh is fixed once for all (for ε = 1, the time computation for the micro-macro model is about two times as large as the CPU time needed for the explicit discretization of the Vlasov equation when the same numerical parameters are considered). If ε is very small, the explicit computation (Vlasov) becomes unreachable while the MM scheme coincides with the LIM model with nongrowing computational cost.

Periodic boundary condition: first order terms in ε

In this subsection, our purpose is to numerically check that the proposed scheme is asymptotically equivalent to the limit scheme, up to O(ε 2 ). To do that, we consider the high-field limit and compare the results given by the numerical scheme of Proposition 3.2 to those obtained by the limit model (3.18). The test case is the bump-on-tail test, the initial condition of which is

f 0 (x, v) = 1 √ 2π 9 10 exp(-v 2 /2) + 2 10 exp(-4(v -4.5) 2 ) (1 + α cos(kx)), where (x, v) ∈ [0, 2π/k]×[-v max , v max ],
with v max = 8, α = 0.05, k = 0.5. The numerical parameters are N x = N v = 256. We plot on Figure 13 the difference between the relative first momentum of g: vg(t, x, v)dv/ε and the corrective terms in the high-field limit model ∂ x ρ -ρE (with N x = 256), as a function of ε. Indeed, according to the previous computations (see subsection 3.2), these two terms only differ from ε 2 terms. We consider the L 1 norm in x to evaluate the difference at the given time t = 3. As expected, we can observe in Figure 13 that the relative difference follows a line of slope 1, which confirms the formal analysis performed in previous sections. We also plot on Figure 14 the two spatial functions vg(t, x, v)dv/ε given by the micro-macro model (MM) and ∂ x ρ(t, x)ρ(t, x)E(t, x) given by the limit model (LIM), at t = 3, for ε = 0.1, 0.08, 0.02, 0.001. It confirms that the scheme (3.2) accurately preserves the high-field limit asymptotics up to O(ε 2 ).

Dirichlet boundary condition and boundary layer

Now, we are interested in the boundary layer induced by an imposed nonequilibrium boundary condition (i.e. which does not correspond to the diffusion limit of the distribution function). Hence, in the whole domain, the solution converges towards the diffusion limit as ε → 0 but when the boundary condition is imposed to be different from the diffusion limit, the solution presents a boundary layer (see [START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | A Numerical Method for Kinetic Semiconductor Equations in the Drift Diffusion Limit[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]). To observe this phenomenon, we consider the following conditions for the kinetic model

f (t = 0, x, v) = 0, f L (t, x = 0, v) = vM (v), f R (t, x = 1, v) = 0, v ∈ R.
For the micro-macro model, we use the procedure presented in section 2.3.3. Finally, we need to determine the consistent boundary condition for the diffusion model; the exact condition is given by [START_REF] Chandrasekhar | Radiative transfer[END_REF] but a good approximation is given by (see [START_REF] Golse | A Numerical Method for Computing Asymptotic States and Outgoing Distributions for a Kinetic Linear Half Space Problem[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | A Numerical Method for Kinetic Semiconductor Equations in the Drift Diffusion Limit[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]) k(t, v) = M (v) w>0 wf L (t, x = 0, w)dw w>0 M (w)dw

Conclusion

The micro-macro decomposition has been applied to collisional Vlasov equations to derive an asymptotic-preserving numerical scheme in the diffusion and high-field limit. The derivation is

Figure 1 :

 1 Figure 1: Micro-macro decomposition in the velocity direction.

(2. 3 )

 3 This formulation is equivalent to the original Vlasov-Poisson-BGK equation (1.1)-(1.2)-(1.5), as precised in the following proposition.

Proposition

. 7 )

 7 This formulation is equivalent to the original Vlasov-Poisson-Fokker-Planck equation (1.1)-(1.2)-(1.6), as precised in the following proposition. Proposition 3.1. (formal

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 23456 Figure 2: Linear Landau damping in the diffusion case: ε = 1. Electric energy as a function of time (up left), and density as a function of space x for t = 1 (up right), t = 2 (bottom left) and t = 10 (bottom right).

Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 789101112 Figure 7: Linear Landau damping in the high-field case: ε = 1. Electric energy as a function of time (up left), and density as a function of space x for t = 1 (up right), t = 2 (bottom left) and t = 10 (bottom right).

εFigure 13 :Figure 14 :

 1314 Figure13: Bump-on-tail test case : L 1 norm of the difference between the relative first momentum of g ( vgdv/ε) and the corrective terms in the high-field limit model ∂ x ρ -ρE. t = 3.

,Figure 15 :Figure 16 :Figure 17 :Figure 18 :Figure 19 :Figure 20 :

 151617181920 Figure 15: Boundary layer test case, example 1: density as a function of the spatial variable at t = 0.4.

  Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Time explicit discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Time implicit discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Derivation of the micro-macro model . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Chapman-Enskog expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Time explicit discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	3 High-field limit
	3.1

1 Introduction 2 Diffusion limit 2.1 Derivation of the micro-macro model . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Chapman-Enskog expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3

small values of ε, (iii) the limit diffusion model (2.5) is also plotted. The comparison is performed on four examples. For the following tests, the time steps are taken as in the Landau case i.e.:

-Vlasov: ∆t = C min εβ x , εβ v , ε 2 , -MM: ∆t = min (Cεβ x / max(0, 1 -Cβ x /ε), Cεβ v / max(0, 1 -Cβ v /ε)),

The CFL number C is chosen equal to C = 0.5.

Example 1

We first investigate the radiative transport equation in which

and the electric field is considered constant null, as in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF]. The initial condition is f (t = 0, x, v) = 0 for the kinetic model. The left boundary condition is f L (v) = v which leads to a boundary layer whereas f R (v) = 0.

On Figure 15, the density is plotted at time t = 0.4 in four configurations: the limit model (with N x = 128), the Vlasov equation (with ε = 10 -2 , N x = 8192 and N v = 64) and the micro-macro model (for ε = 10 -2 and 10 -3 , N x = 256 and N v = 64). As in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF], the four solutions are very close inside the domain. Inside the boundary layer, the micro-macro model gives satisfactoring results.

Example 2

In this example, the velocity domain is [-6, 6] and the equilibrium is

On Figure 16, the density is plotted at time t = 0.1 for different values of ε: 1, 10 -1 and 10 -2 . For the Vlasov model limit model and the micro-macro model, N x = 400 and N v = 32 whereas N x = 128 for the diffusion model. We can observe that the micro-macro model and the reference Vlasov equation give very close results in all the domain. As in the previous case, for small ε, the three solutions are very close inside the domain; the boundary layer is well resolved by the micro-macro model.

Example 3 As in the previous example but here the electric field is considered constant equal to 1, as in [START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF]. The same conclusions as in the previous case are available here.

Example 4

The electric field is considered nonlinearly through the Poisson equation ∂ x E(t, x) = ρ(t, x) = f (t, x, v)dv, with the condition 1 0 E(t, x)dx = -1 as in [START_REF] Duclous | Analysis of a High Order Finite Volume Scheme for the Vlasov-Poisson System[END_REF]. The initial condition is f (t = 0, x, v) = 0 for the kinetic model.

On Figure 18, results obtained by the Vlasov equation and the micro-macro model are plotted with ε = 1 for different times: t = 0.1, 0.2, 1. We can observe that the density given by the two models are very close. On Figure 19, the ε parameter is taken equal to 0.1; the density given by the two models is also nearly superimposed. Figure 20 shows results for ε = 10 -2 : near the right boundary (in which equilibrium boundary condition is imposed), they are nearly superimposed, and at x = 0, a boundary layer can be observed (resolved by the spatial discretization). inspired from previous works in which the methodology has been applied to the radiative transfer equation [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF] or to the BGK equation [START_REF] Benoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF].

However, some specifications are needed here: first, the electric field induces transport in the velocity direction (and consequently an associated discretization), and second, the equilibrium associated to the high-field limit is x-dependent which adds some technical difficulties. Hence, the scheme presented here has the following properties:

• it enjoys the asymptotic-preserving property (up to ε 2 order terms), so that the time discretization can be chosen independent of ε;

• the limiting scheme (for small values of ε) can be either explicit or implicit (using [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF]); obviously, the implicit version enables the use of time steps only induced by the transport;

• boundary layers can be treated.

Some numerical results show the good behaviour of the approach compared to the reference solution in the transport regime (given by the Vlasov equation) and in the diffusion or high-field limit (given by the corresponding asymptotic limit model).

Natural perspective of this work would be its extension to more relevant model from the physical point of view, i.e. using the Landau collision operator. The derivation of an asymptotic-preserving scheme of the Vlasov-Landau equation in the fluid limit seems to be interesting and feasible using this approach and the scheme presented in [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF].