Nicolas Crouseilles

Giacomo Dimarco

Marie-Hélène Vignal

Multiscale schemes for the BGK-Vlasov-Poisson system in the quasi-neutral and fluid limits. Stability analysis and first order schemes. *

Keywords: Collisional Vlasov-Poisson system, quasi-neutral limit, fluiddynamic limit, asymptotic preserving schemes, multiscale, stability analysis

This paper deals with the development and the analysis of asymptotic stable and consistent schemes in the joint quasi-neutral and fluid limits for the collisional Vlasov-Poisson system. In these limits, the classical explicit schemes suffer from time step restrictions due to the small plasma period and Knudsen number. To solve this problem, we propose a new scheme stable for choices of time steps independent from the small scales dynamics and with comparable computational cost with respect to standard explicit schemes. In addition, this scheme reduces automatically to consistent discretizations of the underlying asymptotic systems. In this first work on this subject, we propose a first order in time scheme and we perform a relative linear stability analysis to deal with such problems. The framework we propose permits to extend this approach to high order schemes in the next future. We finally show the capability of the method in dealing with small scales through numerical experiments.

Introduction

The modeling and numerical simulation of plasma phenomena is a very active field of research. Indeed, plasmas underlie numerous technological applications 1 as well as more fundamental researches. Plasma science provides the foundation for present application such as industrial processes (semiconductors, lasers, beam, . . . ) but encompasses also important potential applications such as the generation of energy from fusion. Because of the multiscale character of most plasma phenomena, the numerical simulation of plasmas is still an important challenge for the scientific community. According to the physical context, plasmas can be modeled by means of classes of mathematical models: fluids models and kinetic models. On the one hand, the fluid models make evolve macroscopic quantities such as the density, the temperature or the mean velocity, which depend on time and on the three dimensional space. Fluid models are based on the assumption that the system is close to a thermodynamical equilibrium to be valid [START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF]. On the other hand, kinetic models consider the time evolution of a distribution function which gives the probability of a particle to be in a given state in the six dimensional phase space at a given time [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]. Evidently, for realistic simulations the passage from the three dimensional physical space to the six dimensional one, introduces big challenges for mathematicians, physicists and engineers working on the development of numerical methods. In practice, when dealing with kinetic models, it is necessary to assume hypotheses which reduce the dimensionality of the problem but which increase the difficulty of building numerical schemes and their successive analysis. In any case, there are situations in which kinetic models are mandatory in order to well describe the physical phenomena, while there are situations in which fluid models turn out to be sufficient. In fact, a very typical aspect of these phenomena is the presence of multiple spatial and temporal scales which intervene in different positions and at different times [START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF][START_REF] Radtke | On efficient simulations of multiscale kinetic transport[END_REF]. These behaviors make the construction of numerical methods a real challenge.

In this paper, we propose a new numerical method to solve the kinetic collisional Vlasov-Poisson model of plasma physics. In this model, among the different scales which intervene in the system, there are two important physical length and time scales which should be considered: the Debye length and the electron plasma period. On the one hand, the Debye length measures the typical length at which charge unbalances occur. On the other hand, the electron plasma period gives the typical time oscillations (due to the electrostatic forces) which arise to restore the electric neutrality when charge unbalances occur at the scale of the Debye length [START_REF] Chen | Introduction to plasma physics[END_REF]. When the Debye length and the plasma period are both small compared with macroscopic lengths of interest, the so-called quasi-neutral regime is attained. The plasma appears broadly electrically neutral. Another scale we consider in this paper is the one related to collisions. The magnitude of this phenomenon is commonly measured by the so-called Knudsen number which gives the distance between two successive collisions [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF]. The fluid limit consists in letting the Knudsen number be small compared to the typical scale of the system.

We are interested in studying problems in which the quasi-neutral and the fluid limits may occur simultaneously or not. From the numerical point of view, a classical explicit scheme must resolve the micro scales in order to remain stable and consistent. However, this requires very small time steps and phase space cells. But, on the other side, simulations have to be performed on macroscopic lengths, which makes multiscale models very challenging. Of course, asymptotic models can be derived [START_REF] Jüngel | A hierarchy of hydrodynamic models for plasmas. Quasi-neutral limits in the drift-diffusion equations[END_REF] to describe such regimes, but in situations where both quasi-neutral and non quasi-neutral regimes or both the fluid and the kinetic regimes coexist, different solutions can be find. Hence, domain decomposition approaches [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF][START_REF] Dimarco | An asymptotic preserving automatic domain decomposition method for the Vlasov-Poisson-BGK system with applications to plasmas[END_REF] or hybrid methods [START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equation using particles[END_REF][START_REF] Dimarco | The hybrid moment guided Monte Carlo method for the Boltzmann equation[END_REF][START_REF] Dimarco | Fluid solver independent hybrid methods for multiscale kinetic equations[END_REF][START_REF] Homolle | A low-variance deviational simulation Monte Carlo for the Boltzmann equation[END_REF]] can be adopted. However, the connection of the different models and numerical methods demand specific development as well as the interface identification is not always a simple task to solve. Thus, it seems important to develop numerical methods which can handle multiple regimes simultaneously without the time and space restrictions induced by the small scales. The search for schemes free of such constraints has been the subject of a vast literature. This is precisely the scope for which Asymptotic Preserving (AP) methods have been derived in the past [1, 2, 4, 11-13, 18, 19, 26-29, 40, 45]. These methods are able to overcome these restrictions and automatically degenerate to consistent discretizations of the limiting models when the parameters which characterize the microscopic behaviors goes to zero. However, until now and up to the authors knowledge, AP methods have been developed to treat problems in which only one parameter and so only one microscopic scale and the macroscopic scale were present. Let us also observe that AP methods are the natural candidates for building successively domain decomposition methods [START_REF] Dimarco | An asymptotic preserving automatic domain decomposition method for the Vlasov-Poisson-BGK system with applications to plasmas[END_REF].

The main goals of this work are to develop a new framework for the quasi neutrality problem which permits to construct AP methods which are consistent with this limit, to study the stability of the resulting scheme and then to consider the multiple scale problem characterized by the large variations of the Debye length and Knudsen number. More in details, we present a numerical method which is able to handle at the same time the small Debye length and plasma period and the small Knudsen number for the BGK-Vlasov-Poisson system, allowing stable simulations even when the mesh does not resolve these three scales. In the limit in which these parameters goes to zero the scheme is reduced to a consistent discretization of the corresponding asymptotic model (i.e. the Euler-Poisson system or of the quasi-neutral BGK-Vlasov-Poisson system or again of the quasi-neutral Euler-Poisson system). In other words, the method provides a stable solution in which the plasma oscillations and wave-lengths are filtered out, while classical methods are simply unstable. In this first work, in which we consider the extension of the AP methodologies to the case of multiple scales, we present a first order in time scheme and a relative linear stability analysis which proves stability for small values of the Debye length. The idea is based on two main ingredients: first, the reformulation of the Poisson equation (introduced in [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF] and then used in [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF] and [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF]) and second, the construction of IMEX schemes for collisional kinetic equation (studied for instance in [START_REF] Dimarco | Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations[END_REF][START_REF] Dimarco | Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson system for semiconductors[END_REF]). In additiom, we analyze the state of the art of the schemes which are able to handle the quasi-neutrality constraint and we show that the splitting approach, typical of particle methods as PIC methods [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF], or semi-Lagrangian approaches [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF] may suffer from incompatibility with the quasi-neutrality if no special care is taken. This may suggest that methods which are not based on the splitting approach may be more indicated when we want to increase the order of accuracy. In fact, it seems that the correct use of implicit-explicit schemes for the time discretization of the problem will permit to increase the global order of accuracy of the method here proposed in a near future [START_REF] Crouseilles | Multiscale schemes for the BGK-Vlasov-Poisson system in the quasi-neutral and fluid limits[END_REF].

The paper is organized as follows. In the next Section, we describe the collisional Vlasov-Poisson system in physical and rescaled variables. Then, we recall the reformulated Poisson equation and the quasi-neutral and fluid limits.

In Section 3, we analyze some existing approaches and we introduce our method. Section 4 is devoted to a linear stability analysis while in Section 5 we propose numerical tests which permit to understand the behavior of the scheme. A final Section is devoted to conclusions and to future developments.

2 The BGK-Vlasov-Poisson system and its fluid and quasi-neutral limits

The BGK-Vlasov-Poisson model in physical and rescaled variables

We are interested in the kinetic description of a plasma. For the sake of clarity, we consider a simple one-species model. That means that we assume in the plasma the presence of a uniform static background of ions and we only study the electrons motion described by the Vlasov equation. This equation is coupled with the Poisson equation for the calculation of the electric potential. Thus, the one-species system reads

∂ t f + v • ∇ x f + e m ∇ x ϕ • ∇ v f = Q(f ), ( 1a
) ∆ϕ = e ε 0 (ρ -ρ 0 ), with ρ = f dv, (1b)
where f (x, v, t) is the electron distribution function and where the position and velocity variables x and v are such that (x, v) ∈ Ω × R dv , with Ω ⊂ R dx with d x and d v = 1, 2, or 3 and t > 0 is the time. In the above equations, e is the positive elementary charge, m is the electron mass, ϕ the electric potential, Q(f ) is an operator characterizing the collisions between electrons. Finally, ε 0 is the vacuum permittivity and ρ 0 is the given constant ion density. From now on, we omit the dependence of f on x, v, t unless strictly necessary and we take d x = d v = d even if all the theory naturally extends to the case in which the physical and the velocity spaces have different dimensions. The simplifying choice of considering only one species reduces the system but it permits to maintain its main features, for this reason the extension of the schemes to the general case of multi-species appears to be straightforward.

Let us now precise the collision operator Q(f ). It characterizes the particles interactions and satisfies the following local conservation properties

R d φ(v)Q(f ) dv =: φQ(f ) = 0, (2)
where

φ(v) = 1, v, |v| 2 2 T
are the collision invariants. In this work we use the BGK operator which substitutes the binary interactions of the Boltzmann operator with a relaxation towards the equilibrium. It takes the form (see [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF])

Q BGK (f ) = ν (M [f ] -f ),
where ν = ν(ρ, T ) > 0 is a given relaxation frequency and measures the average time between two collisions while ρ and T are the electron density and temperature defined below. The validity of this operator in describing the physics of non equilibrium phenomena, has been the subject of many papers in the past (see the books [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF] and references therein). In the present paper, we do not discuss it, we just stress that we use this collision operator and that we wish to extend our investigations to more realistic models in the next future. The local Maxwellian equilibrium function takes the form

M [f ] = M (ρ, u, T ) = ρ (2 π k B T /m) d/2 exp -m |v -u| 2 2 k B T ,
where ρ, u, T are the density, mean velocity and temperature of the gas in the x-position and at time t defined as

(ρ, ρu, W ) T = φf = φM [f ] , ρ T = 2 d W - ρ|u| 2 2 . (3)
Now, in order to study the quasi-neutral and fluid limits, we introduce the following scaling of the BGK-Vlasov-Poisson system which highlights the role of the different terms. The scaled variables and unknowns are defined by

x = x x 0 , v = v v 0 , t = v 0 x 0 t, f = (v 0 ) d ρ 0 f, φ = e k B T 0 ϕ, ρ = ρ ρ 0 , ū = u v 0 , T = T T 0 ,
where x 0 > 0 is the typical length of the problem, v 0 = (k B T 0 /m) 1/2 ∈ R is the thermal electron velocity scale with T 0 the electron temperature. We also assume that there exists ν such that ν(ρ, T ) = ν 0 ν(ρ, T ) where ν 0 is the typical relaxation time.

Inserting this scaling into the previous one-species model (1) and omitting the bars, we get the following scaled BGK-Vlasov-Poisson model, denoted by system P ε,λ :

∂ t f + v • ∇ x f + ∇ x ϕ • ∇ v f = ν ε (M [f ] -f ), (4a) λ 2 ∆ϕ = ρ -1, with ρ = f dv, (4b)
with

M [f ] = M (ρ, u, T ) = ρ (2 π T ) d/2 exp -|v -u| 2 2 T , (5)
and the rescaled parameters λ and ε are given by

λ = λ D x 0 = ε 0 k B T 0 x 2 0 e 2 ρ 0 1/2 , ε = v 0 x 0 ν 0 ,
which are called respectively the scaled Debye length and the Knudsen number.

The scaled Debye length is the ratio between the Debye length λ D and the typical length of the problem x 0 . The Knudsen number is the ratio between the typical relaxation time and the typical time of the problem. We recall the Debye length measures the scale of electric interactions in the plasma while the Knudsen number measures the scale of collisions. When the plasma is very dense, both the Debye length and the Knudsen number are very small. However, there are situations in which one scale can be small while the other one is large and vice-versa. In the case of a small Debye length, and so a small plasma period, the plasma appears macroscopically electrically neutral. This is the quasi-neutral regime which can be described by the quasineutral model obtained taking the limit λ → 0 in the scaled BGK-Vlasov-Poisson model P ε,λ . In the case of a small Knudsen number, the plasma appears macroscopically in equilibrium. This is the fluid limit model obtained taking the limit ε → 0 in the scaled BGK-Vlasov-Poisson model P ε,λ . The quasi-neutral limit in plasmas has been theoretically investigated in [START_REF] Brenier | Convergence of the Vlasov-Poisson system to the incompressible Euler equations[END_REF][START_REF] Gérard-Varet | Quasineutral limit of the Euler-Poisson system for ions in a domain with boundaries[END_REF][START_REF] Golse | The Vlasov-Poisson system with strong magnetic field in quasineutral regime[END_REF][START_REF] Kwan | Stability issues in the quasineutral limit of the one-dimensional Vlasov-Poisson equation[END_REF][START_REF] Jüngel | A hierarchy of hydrodynamic models for plasmas. Quasi-neutral limits in the drift-diffusion equations[END_REF] for kinetic or fluid models. The fluid limit has been investigated by different authors in the past, we recall Caflisch and Nishida [START_REF] Caflisch | The fluid dynamical limit of the nonlinear Boltzmann equation[END_REF][START_REF] Nishida | Fluid dynamical limit of the nonlinear Boltzmann equation at the level of the compressible euler equations[END_REF], while for recent theoretical results we refer to the works of Levermore, Golse and Saint-Raymond [START_REF] Golse | The navier-stokes limit of the Boltzmann equation for bounded collision kernels[END_REF][START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF].

Note that in the case of the one-species model, the scaled plasma period equals the scaled Debye length, since

τ = τ p x 0 /v 0 = 1 x 0 /v 0 ε 0 m e 2 ρ 0 1/2 = 1 x 0 ε 0 m v 2 0 e 2 ρ 0 1/2 = λ,
where τ p is the plasma period.

In the next section, we briefly recall the quasi-neutral and fluid limits of the one-species kinetic model (4).

Fluid and quasi-neutral limits of the BGK-Vlasov-Poisson system

Here, we precise the one-species fluid and quasi-neutral models. We start with the fluid limit obtained letting ε tends to 0. Then, we establish the quasi-neutral system with λ → 0, and the joint fluid and quasi-neutral limit corresponding to (ε, λ) → (0, 0). We show that both resulting systems must be reformulated.

2.2.1 Fluid limit: the Euler-Poisson system P 0,λ Formally passing to the limit ε → 0 in (4), we obtain f = M [f ]. Thus, we can show that, at least formally, we recover from the BGK-Vlasov-Poisson system (4), the Euler-Poisson system, denoted by system P 0,λ :

∂ t ρ + ∇ x • (ρ u) = 0, ( 6a
) ∂ t (ρ u) + ∇ x • (ρu ⊗ u) + ∇ x p = ρ ∇ x ϕ, (6b)
∂ t W + ∇ x • ((W + p) u) = ρ u • ∇ x ϕ, (6c)
λ 2 ∆ϕ = ρ -1, (6d)
with

p = 2 d W - ρ |u| 2 2 . ( 7
)
Note that the above limit model can still be obtained for more general collision operators Q(f ) provided they satisfy (2) and admit local equilibrium of the form (5).

2.2.2 Quasi-neutral limit: the quasi-neutral BGK-Vlasov system P ε,0 and its reformulation RP ε,0

Let us turn to the quasi-neutral limit of the BGK-Vlasov-Poisson system. Formally, passing to the limit λ → 0 in (4) merely amounts to replace the Poisson equation (4b) by the quasi-neutrality constraint

ρ = 1. (8)
Then, the quasi-neutral BGK-Vlasov system, denoted by P ε,0 is given by (4a) and [START_REF] Chen | Introduction to plasma physics[END_REF].

In other words, the electrostatic potential becomes the Lagrange multiplier of the constraint [START_REF] Chen | Introduction to plasma physics[END_REF]. This is exactly the same situation as in the incompressible Euler equations in which the pressure is the Lagrange multiplier of the divergence-free constraint. Thus, in this case, in order to recover an explicit equation for the potential ϕ, one idea consists in the reformulation of the system P ε,λ (see [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF] for the fluid case and [START_REF] Degond | An asymptotically stable Particle-In-Cell (PIC) scheme for collisionless plasma simulations near quasineutrality[END_REF], [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF], [START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF] for the kinetic one). Let us begin integrating (4a) with respect to the velocity variable, using the quasi-neutrality constraint, it leads to the divergence-free constraint for the momentum

∇ x • (ρ u) = ∇ x • vf dv = 0. (9)
Now, taking the first moment of the Vlasov equation (4a) yields

∂ t (ρ u) + ∇ x • S = ρ ∇ x ϕ, (10)
where

S = f v ⊗ v dv. (11)
Taking the divergence of equation ( 10) and using the divergence-free constraint obtained from the density equation ( 9), yields an explicit equation for the potential ∇ 2

x : S = ∇ x • (ρ∇ x ϕ). The quasi-neutral BGK-Vlasov system, denoted by system RP ε,0 (for Reformulated P ε,0 ), is thus given by

∂f ∂t + v • ∇ x f + ∇ x ϕ • ∇ v f = ν ε (M [f ] -f ), (12a)
∇ x • (ρ∇ x ϕ) = ∇ 2 x : S. (12b)
Let us observe that the following Lemma holds:

Lemma 2.1 The reformulated BGK-Vlasov-Poisson system RP ε,0 is equivalent to P ε,0 if and only if the initial condition is well prepared to the quasi-neutral regime. More precisely

RP ε,0 =⇒ P ε,0 if and only if ρ(x, 0) = 1, ∇ x • (ρ u)(x, 0) = 0.
Proof: It is sufficient to observe that taking the divergence of the first order moment of equation (12a) we get ∂ t (∇ x • (ρ u)) = 0 and thus

∇ x • (ρ u)(x, t) = ∇ x • (ρ u)(x, 0).
On the other hand, the continuity equation gives

ρ(x, t) = ρ(x, 0) + t ∇ x • (ρ u)(x, 0).
This shows that RP ε,0 gives P ε,0 if and only if ρ(x, 0) = 1 and ∇ x • (ρ u)(x, 0) = 0, which means if and only if the initial condition is well prepared to the quasineutral regime.

2.2.3

Fluid and quasi-neutral limits: the quasi-neutral Euler system P 0,0 and its reformulation RP 0,0

It remains to look at the joint limit (ε, λ) → (0, 0) which can be also formally obtained taking the limit λ → 0 in P 0,λ or ε → 0 in P ε,0 . The limit system P 0,0 is given by (6a)-(6c), [START_REF] Chen | Introduction to plasma physics[END_REF]. Like previously, in this limit we lose the equation for the electric potential and we can recover an equation by applying the same reformulation as before. We get the quasi-neutral Euler system, denoted by

RP 0,0 ∂ t ρ + ∇ x • (ρ u) = 0, (13a) ∂ t (ρ u) + ∇ x • (ρu ⊗ u) + ∇ x p = ρ ∇ x ϕ, ( 13b
) ∂ t W + ∇ x • ((W + p) u) = ρ u • ∇ x ϕ, (13c) ∇ x • (ρ∇ x ϕ) = ∇ 2 x : S, (13d)
with p given by ( 7) and where

S = p Id + ρu ⊗ u.
Like in the previous paragraph, we can prove the following Lemma Lemma 2.2 The reformulated system RP 0,0 is equivalent to the original system P 0,0 if and only if the initial condition is well prepared to the quasi-neutral regime. More precisely

RP 0,0 =⇒ P 0,0 if and only if ρ(x, 0) = 1, ∇ x • (ρ u)(x, 0) = 0.
The proof is identical to that of Lemma 2.1 and we omit it.

2.3 Reformulations of the BGK-Vlasov-Poisson system (RP ε,λ ) and of the Euler-Poisson system (RP 0,λ )

In the previous section, it has been showed that in order to recover an equation for the electric potential in the quasi-neutral limit, one possibility is to reformulate the starting systems either for the BGK-Vlasov-Poisson (P ε,λ ) or for the Euler-Poisson one (P 0,λ ). Thus, in order to derive numerical schemes which work independently on the Debye length scale λ, one idea is to discretize the reformulated system instead of the original one. This is what has been done in [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF], [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Degond | An asymptotically stable Particle-In-Cell (PIC) scheme for collisionless plasma simulations near quasineutrality[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF], in which different numerical schemes efficient in the quasi-neutral limit have been developed. Let us consider the BGK-Vlasov-Poisson system and let us perform the same reformulation as those done for getting the Reformulated BGK-Vlasov system (RP ε,0 ): First integrate (4a) with respect to the velocity variable, this gives the continuity equation. Then, derive with respect to time the continuity equation. Take the divergence of the momentum equation obtained by taking the first order moment of the Vlasov equation (4a). Substract the divergence of the momentum equation to the continuity equation : this yields the reformulated Poisson equation.

Thus, the reformulated BGK-Vlasov-Poisson system, denoted by RP ε,λ , is given by ∂f ∂t

+ v • ∇ x f + ∇ x ϕ • ∇ v f = ν ε (M [f ] -f ), (14a)
λ 2 ∂ tt ∆ϕ + ∇ x • (ρ∇ x ϕ) = ∇ 2 x : S, (14b)
where S is defined by [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF]. In the quasi-neutral limit (λ → 0), the reformulated Poisson equation (14b) formally converges towards the quasi-neutral elliptic equation (12b) satisfied by the potential ϕ. This means that it does not degenerate into an algebraic equation like the Poisson equation (4b) does.

This reformulated system is the appropriate framework to deal with problems which are partly or totally in the quasi-neutral regime. Like previously, the following result can easily be proved.

Lemma 2.3 The reformulated system RP ε,λ gives P ε,λ if and only if the electric potential and its time derivative satisfy the two Poisson equations at the initial time. More precisely

RP ε,λ =⇒ P ε,λ if and only if λ 2 ∆ϕ(x, 0) = ρ(x, 0) -1, λ 2 ∂ t ∆ϕ(x, 0) = -∇ x • (ρ u)(x, 0). ( 15
)
Proof: Taking the divergence of the momentum equation, we get

∂ t (∇ x • (ρ u)) = -∇ 2 x : S + ∇ x • (ρ ∇ x ϕ).
Now using (14b) and integrating with respect to time, yields

∇ x • (ρ u)(x, t) = -λ 2 ∂ t ∆ϕ(x, t) + ∇ x • (ρ u)(x, 0) + λ 2 ∂ t ∆ϕ(x, 0) .
Finally, inserting this result in the continuity equation leads to

ρ(x, t)-1 = λ 2 ∆ϕ(x, t)+ ρ(x, 0) -1 -λ 2 ∆ϕ(x, 0) +t ∇ x • (ρ u)(x, 0) + λ 2 ∂ t ∆ϕ(x, 0) , which proves Lemma 2.3.
The reformulated Euler-Poisson system RP 0,λ which is given by ( 6) and (14b) can be derived in an analogous way. This system is equivalent to P 0,λ if and only if [START_REF] Crouseilles | A forward semi-lagrangian scheme for the numerical solution of the Vlasov equation[END_REF] is satisfied.

Let us conclude this section with two remarks.

Remark 2.4 1-Note that these reformulated systems show that the singularity of the quasi-neutral limit is related to fast oscillations in time at a frequency proportional to 1/λ commonly called the plasma frequency. Indeed, the time behavior of the Laplacian of the electric field is described by the following second order ordinary differential equation with constant coefficients

λ 2 y (t) + ρ y = f.
The exact solution of this equation, is given by y

(t) = A cos √ ρ t/λ +B sin √ ρ t/λ + f /ρ
. This shows that in the quasi-neutral limit, the solution oscillates at a frequency proportional to 1/λ around the quasi-neutral equilibrium. And thus, the singularity of the limit is related to fast oscillations in time.

2-For all λ > 0, equation [START_REF] Crouseilles | A forward semi-lagrangian scheme for the numerical solution of the Vlasov equation[END_REF] gives initial conditions for the reformulated Poisson equation (14b). However, in the limit λ → 0 these initial conditions degenerate into constraints on the initial condition f 0 . These constraints impose that the initial condition f 0 is well prepared to the quasi-neutral regime. This means that if the system RP ε,λ is used, the initial condition must be well prepared to the quasi-neutral regime in the limit. This will be an important remark for the construction of asymptotic preserving schemes.

3 An asymptotic preserving method for the BGK-Vlasov-Poisson system

The construction of numerical schemes capable to capture the fluid and the quasi-neutral limits is closely connected with the notion of asymptotic preserving schemes introduced in [START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF]. In this section, we begin with the definition of asymptotic preserving schemes in the fluid and/or quasi-neutral limit. Then, in section 3.2, we present some classical schemes used for the discretization of the kinetic model and we analyze their asymptotic preserving properties. In section 3.3, we present the state of the art of the AP schemes previously developed in the literature. Then, in sections 3.4, 3.5, we present our new scheme and we prove its uniform stability in the fluid framework.

Definition of the AP properties

In agreement with [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF][START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | An asymptotically stable Particle-In-Cell (PIC) scheme for collisionless plasma simulations near quasineutrality[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF][START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF][START_REF] Pareschi | Implicit-explicit Runge-Kutta methods and applications to hyperbolic systems with relaxation[END_REF], we consider that an asymptotic preserving scheme in a considered limit related to given parameters is a scheme which is both asymptotically stable and consistent.

The asymptotic stability property means that the stability must be proved uniformly with respect to the parameters, that is to say for all choices of the parameters including the limit value 0.

The asymptotic consistency property means that the numerical scheme must be able to approximate all corresponding systems from the original perturbation problem up to the limit systems, in our case RP ε,λ , i.e. (14), RP ε,0 i.e. [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF], RP 0,λ , i.e. (6a)-(6c) and (14b) and RP 0,0 , i.e. [START_REF] Crouseilles | Asymptotic preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF].

Here, we give the following definitions of asymptotic preserving methods for the reformulated system (14) : Definition 3.1 Strong AP property : A consistent and stable discretization method for the system RP ε,λ , i.e. ( 14), of time and space steps ∆t and ∆x, is asymptotic preserving (AP) in the fluid limit if, independently of the initial data, it is stable uniformly with respect to ε for all λ > 0 and if in the limit ε → 0, λ > 0, this discretization becomes a consistent discretization method for the reduced system RP 0,λ , i.e. ( 6)-(14b). It is asymptotic preserving (AP) in the quasi-neutral limit if, independently of the initial data, it is stable uniformly with respect to λ for all ε > 0 and if in the limit λ → 0, ε > 0, this discretization becomes a consistent discretization method for the reduced system RP ε,0 , i.e. [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF]. It is asymptotic preserving (AP) in the joint quasi-neutral and fluid limits if, independently of the initial data, it is stable uniformly with respect to λ and ε and if in the limit (ε, λ) → (0, 0) this discretization becomes a consistent discretization method for the reduced system RP 0,0 , i.e. ( 6)-(12b).

Weak AP property : A weak version of the AP property is a scheme which satisfies the previous strong AP property only if the initial data are consistent with the limit systems.

These properties are summarized in Figure 1.

RP ε,λ RP ε,λ ∆t,∆x ∆t, ∆x → 0 ∆t, ∆x → 0 ε → 0 ε → 0 6 - RP 0,λ RP 0,λ ∆t,∆x 6 - RP ε,λ RP ε,λ ∆t,∆x ∆t, ∆x → 0 ∆t, ∆x → 0 λ → 0 λ → 0 6 - RP ε,0 RP ε,0 ∆t,∆x 6 - RP ε,λ RP ε,λ ∆t,∆x ∆t, ∆x → 0 ∆t, ∆x → 0 (ε, λ) → (0, 0) (ε, λ) → (0, 0) 6 - RP 0,0 RP 0,0 ∆t,∆x 6
-Figure 1: RP ε,λ is the original perturbation problem and RP ε,λ ∆t,∆x its consistent and stable approximation with time and space steps ∆t, ∆x. The AP property means that lim ε→0 RP ε,λ ∆t,∆x = RP 0,λ ∆t,∆x , lim λ→0 RP ε,λ ∆t,∆x = RP ε,0 ∆t,∆x and lim

(ε,λ)→(0,0) RP ε,λ
∆t,∆x = RP 0,0 ∆t,∆x are consistent and stable discretizations of respectively RP 0,λ , RP ε,0 and RP 0,0 . Note that these definitions do not imply that the scheme preserves the order of accuracy in time in the stiff limits ε → 0, λ → 0 or (ε, λ) → (0, 0). In this latter case we say that the scheme is asymptotically accurate (AA).

The development of AP and AA schemes for kinetic models in the fluid limit has been already successfully treated in [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Pieraccini | Implicit-Explicit schemes for BGK kinetic equations[END_REF] for the BGK equation and in [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF][START_REF] Dimarco | Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations[END_REF] for the full Boltzmann operator.

Concerning, the asymptotic properties in the quasi-neutral limit as shown later in the Section even if some good and very interesting attempts have been recently done [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF][START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF][START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF], they still present some drawbacks and for this reason a new method which try to overcome these drawbacks will be proposed in Section 3.4.

Classical existing schemes for the quasi-neutral kinetic model

In this section, we analyze the consistency of different classical schemes for the quasi-neutral limit system [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF]. We consider the classical first order splitting and the second order Strang splitting which are often used together with PIC schemes [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF][START_REF] Langdon | Direct implicit large time-step particle simulation of plasmas[END_REF] or semi-Lagrangian methods [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Crouseilles | Conservative semilagrangian schemes for Vlasov-type equations[END_REF][START_REF] Crouseilles | A forward semi-lagrangian scheme for the numerical solution of the Vlasov equation[END_REF]. In this analysis, we omit the relaxation term, which means we consider only the non collisional case, the collisional case gives the same results. In fact, the relaxation operator, conserving the moments of the solution, does not play any role in this situation.

Initially, we start with a quasi-neutral state and we want to see if the schemes are able to propagate the quasi-neutrality or if they introduce a consistency error on this quasi-neutral state.

First order splitting

It is well known that if the first order splitting method is employed for solving only the kinetic equation without coupling with the Poisson equation, the chronological order of appearance of the operators (transport, relaxation, force) does not change the consistency of the scheme. However, this is not the case when the kinetic equation is coupled with another equation like in the Vlasov-Poisson system.

The initial conditions satisfy the quasi-neutral limit, i.e. ρ n = 1 and ∇ x • (ρu) n = 0. We solve first the transport part and then the force part for system [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF] with ν = 0. We discretize only the time variable for the distribution f , we keep the phase space continuous and we do not put intentionally any superscript to the electric potential for the moment, we discuss after the best choice. We have

f * = f n -∆t v • ∇ x f n , f n+1 = f * -∆t ∇ x ϕ • ∇ v f * .
Then, integrating over the velocity space, we get ρ * = ρ n -∆t∇ x • (ρu) n = 1, (ρu) * = (ρu) n -∆t∇ x • S n , from which we deduce ∇ x • (ρu) * = -∆t∇ 2

x : S n .

The second step of the splitting gives for the density and the momentum ρ n+1 = ρ * = 1, (ρu) n+1 = (ρu) * + ∆tρ * ∇ x ϕ, from which we deduce the divergence of the momentum

∇ x • (ρu) n+1 = ∇ x • (ρu) * + ∆t∇ x • (ρ * ∇ x ϕ) = ∆t(-∇ 2 x : S n + ∇ x • (ρ n ∇ x ϕ))
. Then, we ensure that ∇ x • (ρu) n+1 = 0 and ρ n+1 = 1 if we choose the electric potential satisfying the quasi-neutral equation ∇ 2

x :

S n = -∇ x • (ρ n ∇ x ϕ).
And, the quasi-neutral state is preserved for all times. Now, we analyze the same first order splitting but, in which we solve first the force term and then the transport term. The initial condition satisfies again the quasi-neutral conditions, i.e. ρ n = 1 and ∇ x • (ρu) n = 0. We have

f * = f n -∆t∇ x ϕ • ∇ v f n , f n+1 = f * -∆t v • ∇ x f * .
Integrating over the velocity space gives for the density and the momentum for the intermediate step

ρ * = ρ n = 1, (ρu) * = (ρu) n + ∆tρ n ∇ x ϕ so that ∇ x • (ρu) * = ∆t∇ x • (ρ n ∇ x ϕ). The second step gives ρ n+1 = ρ * -∆t∇ x • (ρu) * = ρ n -∆t 2 ∇ x • (ρ n ∇ x ϕ), (ρu) n+1 = (ρu) * -∆t∇ x • S * . And, then ∇ x • (ρu) n+1 = ∇ x • (ρu) * -∆t∇ 2 x : S * = ∆t(∇ x • (ρ n ∇ x ϕ) -∇ 2
x : S * ). This means that in principle, we can choose an electric potential which ensures that

∇ x • (ρu) n+1 = 0, this is the case if ϕ is the solution of ∇ x • (ρ n ∇ x ϕ) = ∇ 2
x : S * . However, in the general case, there is no choice of the electric potential which ensures the propagation of the quasi-neutral state ρ n+1 = ρ n = 1.

Let us observe that both methods proposed in [START_REF] Degond | An asymptotically stable Particle-In-Cell (PIC) scheme for collisionless plasma simulations near quasineutrality[END_REF][START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF] and [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF] are based on the same time splitting for the solution of the Vlasov equation and then on a discretization of the phase space by, in the first case, a particle approach while in the second case by means of semi-Lagrangian techniques. These approaches introduce a consistency error at the limit λ → 0 of order ∆t 2 at each time step.

Second order Strang splitting

We analyze here the second order Strang splitting. From the previous analysis, it appears that with an order one splitting, it is sufficient to make a good choice of the chronological order of each operators to preserve the quasi-neutrality. But here, we will see that when a second order splitting is applied, there is no choice of the first operator which preserves the propagation of the quasi-neutral states.

We start again with an initial condition which satisfies the quasi-neutral limit, i.e. ρ n = 1 and ∇ x • (ρu) n = 0. We begin solving the transport part for system [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF] with ν = 0. We have

f * = f n - ∆t 2 v • ∇ x f n , f * * = f * -∆t∇ x ϕ • ∇ v f * , f n+1 = f * * - ∆t 2 v • ∇ x f * * .
Then, integrating over the velocity space we get ρ * = ρ n -∆t 2 ∇ x • (ρu) n = 1, (ρu) * = (ρu) n -∆t 2 ∇ x • S n , and for the second step we get ρ * * = ρ * = 1, (ρu) * * = (ρu) * + ∆tρ * ∇ x ϕ, from which we obtain, taking the divergence of the momentum

∇ x •(ρu) * * = ∇ x •(ρu) * +∆t∇ x •(ρ * ∇ x ϕ) = ∆t(-∇ 2 x : S n +∇ x •(ρ n ∇ x ϕ)). (16)
Finally, the density for the third step becomes

ρ n+1 = ρ * * - ∆t 2 ∇ x • (ρu) * * , (ρu) n+1 = (ρu) * * + ∆t 2 ∆tρ * * ∇ x ϕ.
Now, we can consider an electric potential such that ∇ x • (ρ n ∇ x ϕ) = ∇ 2 x : S n , thanks to equation ( 16), we have ∇ x • (ρu) * * = 0 and then ρ n+1 = 1. However, there is no choice which ensures ∇ x • (ρu) n+1 = 0. This implies that at the next time step we lose the quasi-neutral state. If now, we repeat the above computation starting from the force term instead of the transport term in the Strang splitting, we end with the conclusion ρ n+1 = 1 and then quasineutrality is again lost. Thus, it appears that there is no easy solution for the construction of high order schemes preserving quasi-neutral states for system [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF] when splitting methods are used.

State of the art for the fluid model in the quasi-neutral limit

AP schemes in the quasi-neutral limit for the Euler-Poisson system have been proposed in [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF], [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF] and [START_REF] Vignal | A boundary layer problem for an asymptotic preserving scheme in the quasi-neutral limit for the Euler-Poisson system[END_REF]. All these schemes are based on a reformulation of the Poisson equation and on an implicit treatment of the mass flux term

∇ x • (ρ u).
We present the idea in the simple case of the isentropic one-fluid-Euler-Poisson system:

∂ t ρ + ∇ x • (ρ u) = 0, ∂ t (ρ u) + ∇ x • S = ρ ∇ x ϕ, λ 2 ∆ϕ = ρ -1, (17)
where

S = ρ u ⊗ u + C ρ γ .
Then, the AP scheme proposed in [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF], [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF] and [START_REF] Vignal | A boundary layer problem for an asymptotic preserving scheme in the quasi-neutral limit for the Euler-Poisson system[END_REF] and studied in [START_REF] Degond | Analysis of an asymptotic preserving scheme for the euler-poisson system in the quasineutral limit[END_REF], is the following

ρ n+1 -ρ n ∆t + ∇ x • (ρ u) n+1 = 0, ( 18a
) (ρ u) n+1 -(ρ u) n ∆t + ∇ x • S n = ρ n ∇ x ϕ n+1 , (18b)
λ 2 ∆ϕ n+1 = ρ n+1 -1. (18c)
In these works, it is proved that the above scheme has an uncoupled formulation. Indeed, the discretized Poisson equation can be rewritten using density and momentum equation as

λ 2 ∆ϕ n+1 -λ 2 ∆ϕ n ∆t - λ 2 ∆ϕ n -λ 2 ∆ϕ n-1 ∆t = ρ n+1 -ρ n ∆t - ρ n -ρ n-1 ∆t = -∇ x • (ρ u) n+1 -∇ x • (ρ u) n = ∆t ∇ 2 x : S n -∆t ∇ x • (ρ n ∇ x ϕ n+1 ).
Thus, for all n ≥ 2 we get

λ 2 ∆ϕ n+1 -λ 2 ∆ϕ n ∆t - λ 2 ∆ϕ n -λ 2 ∆ϕ n-1 ∆t ∆t + ∇ x • (ρ n ∇ x ϕ n+1 ) = ∇ 2 x : S n , (19)
where ϕ 0 and ϕ 1 can be computed by solving the discrete Poisson equation λ 2 ∆ϕ n+1 = ρ n+1 -1. Let us observe that this equation is nothing else but a discretization of the reformulated Poisson equation (14b). This means that in the limit this scheme is consistent with the reformulated Euler-Poisson system RP 0,0 . Moreover, we can easily check that if the initial data are consistent with the quasi-neutral limit, the scheme remains consistent for all t > 0.

In [START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF] it is proposed an alternative to the above scheme. This alternative permits to get a scheme consistent with the quasi-neutral limit even if the initial data are not prepared to this regime. It consists in inserting successively the mass and momentum equations in the Poisson equation, to get

λ 2 ∆ϕ n+1 = ρ n+1 -1 = ρ n -1 -∆t ∇ x • (ρ u) n+1 = ρ n -1 -∆t ∇ x • (ρ u) n + (∆t) 2 ∇ 2 x : S n -∆t 2 ∇ x • (ρ n ∇ x ϕ n+1 ),
which can be rewritten as

∇ x • (λ 2 + ∆t 2 ρ n ) ∇ x ϕ n+1 = ρ n -1 -∆t ∇ x • (ρ u) n + ∆t 2 ∇ 2 x : S n . ( 20
)
This last relation gives an uncoupled scheme which does not need any initial resolution of the discrete Poisson equation. This last reformulation has been proposed in [START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF] in the context of the Vlasov equation. However, while it is used for solving the simplified Euler-Poisson system (17), it gives a scheme which gives exactly the quasi-neutral constraint in the limit λ → 0 independently from the initial data. This is not the case for the Vlasov equation. In fact, as already explained, in this latter case, the splitting between the transport and the force terms introduces a consistency error on the quasi-neutral constraint in the limit λ → 0. Moreover, in the fluid case, the implicit treatment of the mass flux term is no more expensive than a classical explicit scheme, while if we want to bring the scheme [START_REF] Degond | An asymptotically stable Particle-In-Cell (PIC) scheme for collisionless plasma simulations near quasineutrality[END_REF] with the modified discretized Poisson equation [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF] to the kinetic case this requires the resolution of large linear systems which in realistic situations should be avoided. For these reasons, in the next section, we present an AP scheme in which the flux terms are treated explicitly and which is consistent with the quasi-neutral limit.

A new asymptotic preserving scheme in the quasineutral limit

In this section, we introduce the new numerical scheme for the BGK-Vlasov-Poisson system. We have seen that for the Euler-Poisson case the schemes proposed in [START_REF] Crispel | An asymptotically stable discretization for the Euler-Poisson system in the quasineutral limit[END_REF], [START_REF] Crispel | An asymptotically preserving scheme for the two-fluid Euler-Poisson model in the quasineutral limit[END_REF] seem to be the good choice for quasi neutrality problems.

However, If we want to extend this approach to the Vlasov equation, we must apply an implicit treatment to the space transport term in the Vlasov equation.

Then, the collisionless Vlasov equation should be discretized as follows

f n+1 -f n ∆t + v • ∇ x f n+1 + ∇ x ϕ n+1 • ∇ v f n = 0,
where the electric potential can be computed by either the discrete reformulated Poisson equation ( 19) or [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF]. In this way, taking the velocity moments, we obtain a discretization of the continuity equation with an implicit mass flux term but also a momentum equation with an implicit flux term. Actually, there is no discretization of the Vlasov equation which leads to an implicit treatment of the mass equation and an explicit treatment of the momentum equation. Then, the equivalence of the discrete Poisson equation and of the discrete reformulated Poisson equation is not ensured. Furthermore, as already explained, the implicit treatment of the transport term induces the resolution of a linear system for each discrete velocity and turns in an enormous computational cost for each iteration which should be avoided.

In the present paper, we consider an alternative approach. Let us first consider the quasi-neutral system RP ε,0 , i.e. [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF]. Then, we assume quasineutral initial conditions i.e. ρ n = 1 and ∇ x • (ρu) n = 0, and we write the following time semi-discretization

f n+1 = f n -∆t v • ∇ x f n -∆t∇ x ϕ n+1 • ∇ v f n + ∆t ν ε (M [f n+1 ] -f n+1 ), (21a) ∇ x • (ρ n ∇ x ϕ n+1 ) = ∇ 2 x : S n . ( 21b
)
Note that the implicit treatment of the relaxation source term is the basis of AP schemes in the fluid limit developed in [START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF][START_REF] Dimarco | Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Pieraccini | Implicit-Explicit schemes for BGK kinetic equations[END_REF]. Taking the velocity moments of (21a) leads to

ρ n+1 = ρ n -∆t∇ x • (ρu) n , ( 22a
) (ρu) n+1 = (ρu) n -∆t∇ x • S n + ∆tρ n ∇ x ϕ n+1 , (22b)
which gives

ρ n+1 = ρ n , and ∇ x • (ρu) n+1 = -∆t(∇ 2 x : S n -∇ x • (ρ n ∇ x ϕ n+1 )) = 0.
This proves that the quasi-neutrality constraint is propagated in time. Concerning the fluid limit, solving for f n+1 in the limit ε → 0, one gets

f n+1 = M [f n+1 ],
which implies that the method is AP also in the fluid limit. Let us observe that for computing the Maxwellian state at time n + 1, we need the knowledge of the moments at time n + 1, and the moments of the distribution function are computed explicitly. Extension of this approach to high order schemes for system RP ε,0 by using Implicit-Explicit Runge-Kutta schemes (IMEX) will be discussed in future works.

Now, let us turn to the reformulated BGK-Vlasov-Poisson system RP ε,λ , i.e. [START_REF] Crouseilles | Conservative semilagrangian schemes for Vlasov-type equations[END_REF]. We propose the following time semi-discretization given by (21a) and

λ 2 ∆ϕ n+1 -2∆ϕ n + ∆ϕ n-1 ∆t + ∆t ∇ x • ρ n ∇ x ϕ n+1 = ∆t ∇ 2 x : S n , ∀ n ≥ 2. ( 23
)
Let us observe that the resolution of the electric potential equation requires the knowledge of the potential field at time n and n -1. Then, initially two resolutions of the constrained Poisson equation must be done. In [START_REF] Degond | Asymptoticpreserving particle-in-cell method for the Vlasov-Poisson system near quasineutrality[END_REF], the authors explain that these two initial steps can introduce instabilities in the numerical results.

To bypass this limitation, let us first remark that using the moment equations ( 22), the discrete reformulated Poisson equation ( 23) can be rewritten

λ 2 ∆ϕ n+1 -2∆ϕ n + ∆ϕ n-1 ∆t - ρ n+2 -2 ρ n+1 + ρ n ∆t = 0.
So, forgetting the initial steps, equation ( 23) is equivalent to the following discrete Poisson equation:

λ 2 ∆ϕ n+1 = ρ n+2 -1. ( 24
)
But, now using the moment equations [START_REF] Dimarco | The hybrid moment guided Monte Carlo method for the Boltzmann equation[END_REF], at time n + 2 for (22a) and n + 1 for (22b), this discrete Poisson equation is equivalent to

λ 2 ∆ϕ n+1 = ρ n+1 -1 -∆t∇ x • (ρu) n+1 , = ρ n+1 -1 -∆t∇ x • (ρu) n + ∆t 2 ∇ 2 x : S n -∆t 2 ∇ x • ρ n ∇ x ϕ n+1 .
This gives an AP discretization of the BGK-Vlasov-Poisson system given by

f n+1 = f n -∆t v • ∇ x f n -∆t∇ x ϕ n+1 • ∇ v f n + ∆t ν ε (M [f n+1 ] -f n+1 ), (25a) ∇ x • (λ 2 + ∆t 2 ρ n )∇ x ϕ n+1 = ρ n+1 -1 -∆t∇ x • (ρu) n + ∆t 2 ∇ 2
x : S n . (25b) Let us now analyze the behavior of the above scheme. First, if we let λ → 0 in (25b) we immediately get, if the initial data are consistent with the quasi neutral limit, equation (21b) which means that the scheme is AP in the quasi neutral limit. On the other hand if λ = 0 at t = 0, but the initial data are not consistent with the quasi neutral limit we get after the first time step

ρ 1 = ρ 0 -∆t∇ x • (ρu) 0 , (ρu) 1 = (ρu) 0 -∆t∇ x • (S) 0 + ∆tρ 0 ∇ x ϕ 1 ,
and

∇ x • (ρu) 1 = ∇ x • (ρu) 0 -∆t∇ 2 x : S 0 + ∆t ∇ 2 x : S 0 + 2ρ 1 -ρ 0 -1 ∆t 2 , = ∇ x • (ρu) 0 -∇ x • (ρu) 0 + ρ 1 -1 ∆t = ρ 1 -1 ∆t .
which in the second time step leads to

ρ 2 = ρ 1 -∆t∇ x • (ρu) 1 = ρ 1 -∆t ρ 1 -1 ∆t = 1, (ρu) 2 = (ρu) 1 -∆t∇ x • (S) 1 + ∆tρ 0 ∇ x ϕ 2 ,
which leads to

∇ x • (ρu) 2 = ∇ x • (ρu) 1 -∆t∇ 2 x : S 1 + ∆t ∇ 2 x : S 1 + 2ρ 2 -ρ 1 -1 ∆t 2 , = ∇ x • (ρu) 1 -∇ x • (ρu) 1 + ρ 2 -1 ∆t = ρ 2 -1 ∆t = 0,
and finally to ρ n = 1, ∀n ≥ 2 and ρ n = 1, ∀n ≥ 2. Moreover, if at a given instant of time t n , the Debye length becomes zero, the above analysis can be repeated which shows that at the time step t n+2 the quasi neutrality is obtained and then propagated for all times. Thus, the scheme gives the demanded AP property for the quasi neutral case for both consistent and not consistent initial data. The consistency with the fluid limit remains unchanged and it follows from the fact that f is projected on the relative equilibrium state M [f ] when ε → 0. To conclude this section we prove a linear stability result for above described scheme in the fluid limit ε → 0. This analysis prove that the scheme proposed is stable for all values of λ to small perturbations of the quasi neutral equilibrium state. The algorithm can be summarized as follows, from an initial condition f 0

• advance (22a) to get ρ n+1 and use the same time scheme to get (ρu) n+1 , W n+1 ,

• compute ϕ n+1 from (25b),

• compute (ρu) n+1 , W n+1 ,

• compute M [f n+1 ] from (5)
• advance (25a) to get f n+1 .

3.5 Linear stability analysis of scheme (25) in the fluid limit

The proposed scheme [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] for the one dimensional isentropic Euler-Poisson system reads

ρ n+1 -ρ n ∆t + ∂ x (ρ u) n = 0, (ρ u) n+1 -(ρ u) n ∆t + ∂ x S n = ρ n ∂ x ϕ n+1 , λ 2 ∂ 2 xx ϕ n+1 -2 ∂ 2 xx ϕ n + ∂ 2 xx ϕ n-1 ∆t + ∆t ∂ x ρ n ∂ x ϕ n+1 = ∆t ∂ 2 xx S n , ∀ n ≥ 1,
with S = ρu ⊗ u + p and p = ρ γ , γ > 1 the pressure. Now, linearizing this system around the steady state ρ = 1, q = ρ u = 0 and ∂ x ϕ = 0, we get

               ρ n+1 -ρ n ∆t + ∂ x q n = 0, q n+1 -q n ∆t + ∂ x ρ n = ∂ x ϕ n+1 , λ 2 ∂ 2 xx ϕ n+1 -2 ∂ 2 xx ϕ n -∂ 2 xx ϕ n-1 ∆t 2 + ∂ 2 xx ϕ n+1 = ∂ 2 xx ρ n .
We then apply the spatial Fourier transform on the variables ρ, q and ϕ and we denote these transformed variables by ρ, q and φ. Setting

ψ n+1 = (ϕ n+1 - ϕ n )/∆t, we get                ρn+1 -ρn + i k ∆t qn = 0, qn+1 -qn + i k ∆t ρn -i k ∆t φn+1 = 0, φn+1 -φn -∆t ψ n+1 = 0, λ 2 ψn+1 -λ 2 ψn + ∆t φn+1 -∆t ρn = 0.
In [START_REF] Degond | Analysis of an asymptotic preserving scheme for the euler-poisson system in the quasineutral limit[END_REF], it has been remarked that a system of this kind, i.e. semi-discretized in time, from the stability point of view gives analogous results of a full discretization where central numerical derivatives are employed for the space derivatives. Now, it is well known that this choice leads to an unstable scheme, thus to overcome this problem, following [START_REF] Degond | Analysis of an asymptotic preserving scheme for the euler-poisson system in the quasineutral limit[END_REF], we introduce some numerical viscosity on the fluid equations which depends on a constant α, which leads to the following modified scheme

               ρn+1 -ρn + i k ∆t qn + α k 2 ∆t ρn = 0, qn+1 -qn + i k ∆t ρn + α k 2 ∆t qn -i k ∆t φn+1 = 0, φn+1 -φn -∆t ψ n+1 = 0, λ 2 ψn+1 -λ 2 ψn + ∆t φn+1 -∆t ρn = 0,
which can be rewritten as

B     ρn+1 qn+1 φn+1 ψn+1     = A     ρn qn φn ψn     , with B =     1 0 0 0 0 1 -i k ∆t 0 0 0 1 -∆t 0 0 ∆t λ 2     and A =     1 -α k 2 ∆t -i k ∆t 0 0 -i k ∆t 1 -α k 2 ∆t 0 0 0 0 1 0 ∆t 0 0 λ 2     .
Now, in order to prove stability, we want to prove that the modulus of the eigenvalues µ of B -1 A are lower than 1. Thus, observing that det(B) = λ 2 + ∆t 2 = 0, the dispersion relation (or the characteristic polynomial) reads as

µ 4 -2 λ 2 λ 2 + ∆t 2 + 1 -α k 2 ∆t µ 3 + + λ 2 λ 2 + ∆t 2 + (1 -α k 2 ∆t) 2 + 4 (1 -α k 2 ∆t) λ 2 λ 2 + ∆t 2 + k 2 ∆t 2 λ 2 λ 2 + ∆t 2 µ 2 - 2 λ 2 λ 2 + ∆t 2 1 -α k 2 ∆t + (1 -α k 2 ∆t) 2 + k 2 ∆t 2 µ+ + λ 2 λ 2 + ∆t 2 (1 -α k 2 ∆t) 2 + k 2 ∆t 2 = 0
Let us observe that a necessary condition for the stability is that the constant term (product of the eigenvalues) is in modulus smaller than 1. This condition is ensured, for instance, taking

k 2 ∆t 2 ≤ 1, and (1 -α k 2 ∆t) 2 = 0.
The first condition corresponds to the classical CFL condition of an explicit scheme, i.e. ∆t ≤ 1/k, because 1/k plays the role of the space mesh size ∆x. The second condition corresponds to take a viscosity of the order of α ∼ 1/k ∼ ∆x if the CFL condition is satisfied.

If we consider these conditions fulfilled, the dispersion relation becomes

µ 4 -2 λ 2 λ 2 + ∆t 2 µ 3 + λ 2 λ 2 + ∆t 2 1 + k 2 ∆t 2 µ 2 + -2 λ 2 λ 2 + ∆t 2 k 2 ∆t 2 µ + λ 2 λ 2 + ∆t 2 k 2 ∆t 2 = 0.
We denote then by R and C the following quantities

R = λ 2 λ 2 + ∆t 2 ∈ [0, 1], and C 2 = k 2 ∆t 2 ∈ [0, 1].
Let us now suppose that ∆t > λ, which means that the time steps are larger than the Debye length which is exactly what we ask to our scheme. In this case, R lies in R ∈ [0, 1/2[. We are now able to determine the roots of the dispersion relation by using the Descartes method.

The dispersion relation can be put in the form

P 1 (µ) = µ 4 -2 R µ 3 + R (1 + C 2 ) µ 2 -2 R C 2 µ + R C 2 = 0.
Setting z = µ -R/2, the above polynomial can be rewritten

P 2 (z) = z 4 + p z 2 + q z + r = 0, with p = λ 2 λ 2 + ∆t 2 - 3 2 λ 2 λ 2 + ∆t 2 + 1 + k 2 ∆t 2 = R - 3 2 R + 1 + C 2 , q = λ 2 λ 2 + ∆t 2 - λ 2 λ 2 + ∆t 2 2 + λ 2 λ 2 + ∆t 2 1 + k 2 ∆t 2 -2 k 2 ∆t 2 = = R -R 2 + R 1 + C 2 -2 C 2 , r = λ 2 λ 2 + ∆t 2 - 3 16 λ 2 λ 2 + ∆t 2 3 + 1 4 λ 2 λ 2 + ∆t 2 2 (1 + k 2 ∆t 2 ) - λ 2 λ 2 + ∆t 2 k 2 ∆t 2 + k 2 ∆t 2 = R - 3 16 R 3 + 1 4 R 2 (1 + C 2 ) -R C 2 + C 2 .
Now, we search for real coefficients a, b and c such that

z 4 +p z 2 +q z +r = (z 2 +a z +b) (z 2 -a z +c) = z 4 +(b+c-a 2 ) z 2 +a (c-b) z +b c. (26)
In order to do that, we must solve the non linear system

   b + c = p + a 2 a (c -b) = q b c = r ⇔    2 a b = a (p + a 2 ) -q 2 a c = a (p + a 2 ) + q 2 a b × 2 a c = 4 a 2 r ⇔            b = 1 2 p + a 2 - q a c = 1 2 p + a 2 + q a a 2 (p + a 2 ) 2 -q 2 = 4 a 2 r
Denoting by y = a 2 , the last equation can be recast in the following third order polynomial equation

P 3 (y) = y 3 + 2 p y 2 + (p 2 -4 r) y -q 2 = 0. ( 27
)
One can remark that the above polynomial is negative at 0 and its limit is +∞ at +∞. This means that P 3 admits a real and positive root y 0 and we get a real a = √ y 0 .

Then, we use the Cardan formulas to solve equation [START_REF] Dimarco | Numerical methods for kinetic equations[END_REF]. We start by defining y = ξ -2 p/3 to get

P 4 (ξ) = ξ 3 + p ξ + q = 0, with p = - p 2 3 -4 r = - R 2 3 (1 -C 2 ) 2 -4 R C 2 1 - 2 3 R ≤ 0,
and

q = - 2 p 3 3 3 + 8 r p 3 -q 2 = - 2 27 R 3 (1 + C 2 ) 3 + 4 3 C 2 (2 -C 2 ) R 2 (1 -R).
The discriminant of P 4 (ξ) is given, after some algebra, by

∆ = -4 (p ) 3 27 -(q ) 2 = 16 27 R 3 C 2 (1-R) R 2 (1+11 C 2 -C 4 )+R (C 6 -12 C 4 -8 C 2 )+16 C 4 .
Finally, the polynomial of degree two in the R variable admits for discriminant

∆ = (C 6 -12 C 4 -8 C 2 ) 2 -4 3 C 4 (1 + 11 C 2 -C 4 ) = C 6 (C 6 -24 C 4 + 192 C 2 -512) = C 6 (C 2 -8) 3 < 0 for C ∈ [0, 1]
. This implies ∆ > 0 and that P 4 has 3 real roots. These roots are given by

ξ k = j k -q + i √ ∆ 2 1/3 + j -k -q -i √ ∆ 2 1/3
, for k = 0, 1, 2, with j = exp (2 iπ/3). Then

ξ 0 = 2 Re   -q + i √ ∆ 2 1/3   , and y 0 = ξ 0 - 2 p 3 .
An easy computation shows that y 0 = 0 if R = 0 and on Figure 2, we verify numerically that y 0 is positive for all C ∈ [0, 1] and R ∈ [0, 1/2]. Then, we have a = √ y 0 and from system (3.5)

b = 1 2 p + (a) 2 - q a , and c = 1 2 p + (a) 2 + q a .
Coming back to the dispersion relation we get thanks to (26) Thus, finally we get the following sufficient conditions for stability to be guarantee

P 1 (µ) = P 2 (z) = (z 2 + a z + b) (z 2 + a z + c) = µ 2 + (a -R) µ + b + R 4 -a R 2 µ 2 + (a -R) µ + c + R 4 -a R 2 .
β 1 = b + R 4 -a R 2 ∈ [0, 1],
and

β 2 = c + R 4 -a R 2 ∈ [0, 1], ∆ 1 = (a -R) 2 -4 b + R 4 -a R 2 ≤ 0, and ∆ 2 = (a -R) 2 -4 c + R 4 -a R 2 ≤ 0,
which we verify numerically in Figures 3 and4.

Then, we obtain that the only necessary and sufficient conditions for stability to be ensured are k 2 ∆t 2 ≤ 1, and (1 -α k 2 ∆t) 2 = 0.

We recall that the first condition corresponds to the CFL condition of an explicit scheme, and the second condition corresponds to the choice of the viscosity.

Numerical tests

In this Section, we present one dimensional numerical results for the kinetic model in the fluid and quasi-neutral limits. In particular, we focus our attention on the behavior of this new algorithm in the quasi-neutral limit. We compare our scheme [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] to the classical scheme for the quasi neutrality given by

f n+1 = f n -∆t v ∂ x f n -∆t ∂ x ϕ n+1 ∂ v f n + ν ε (M [f n+1 ] -f n+1 ), ( 28a
)
λ 2 ∂ 2 xx ϕ n+1 = ρ n+1 -1, (28b)
the rest of the scheme, in particular the part concerning collisions, being equal to [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF]. For all test cases, we consider a uniform Eulerian discretization in space and velocity. We denote by ∆x and ∆v the constant space and velocity steps. The velocity space is truncated and replaced by [-v

max -∆v/2, v max + ∆v/2] = ∪ Nv k=1 [v k -∆v/2, v k + ∆v/2]
We use staggered grids for the distribution function and its moments and the electric potential. Then [0, L[= ∪ Nx i=1 [x i -∆x/2, x i + ∆x/2[ is the space mesh; the quantities (f , ρ, u and S) are located at the center of these cells whereas the electric potential is located at the interfaces

x i-1/2 = x i -∆x/2.
Classical Lax-Friedrichs discretizations for each transport operators (in space and in velocity) are employed. The corresponding CFL condition which is the only quantity which limits the time step in our scheme is

∆t ≤ 1 v max ∆x + max i |E n+1 i | ∆v , ( 29
)
where

E n+1 i = -(ϕ n+1 i+1/2 -ϕ n+1 i-1/2
)/∆x is an approximation of the electric field E(x, t) = -∂ x ϕ(x, t) for x in the cell M i = [x i -∆x/2, x i + ∆x/2[ and for t ∈ (t n , t n + ∆t]. Since this condition is non linear, we choose ∆t at time n using the electric field at time n, we update the equations with this time step and we iterate if the condition is not verified with a smaller time step.

We present two different test cases. Both consists in a periodic perturbation of stationary solution of the system. The results prove that our new scheme [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF] is AP in the joint fluid and quasi-neutral limit while the classical scheme (28) is not.

Periodic perturbation of a stationary solution

In this test case, we perturb the moments of a Maxwellian stationary solution. The domain in space is [0, 1] and the initial condition is given by

f (x, v, 0) = ρ 0 (x) √ 2 π exp - (v -u 0 (x)) 2 2 ,
where

ρ 0 (x) = 1 + δ ρ cos(2 π x), ρ 0 (x) u 0 (x) = δ q cos(2 π x).
with δ q = 0.05 and δ ρ = λ/2. Then, the initial condition is not well-prepared to the quasi-neutral regime. For all simulations, N v = N x = 128 and v max = 6 for λ = 1, and v max 18 for λ = 10 -5 . We consider periodic boundary conditions in space and null boundary conditions in velocity. Furthermore, we set homogenous boundary condition for the electric potential. On Figures 5,6, 7, we compare the results given by the classical and AP schemes [START_REF] Dimarco | Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson system for semiconductors[END_REF] and [START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF].

On Figures 5 and6, we can see that the schemes give same results for large Debye-lengths (λ = 1) and small or large Knudsen numbers (ε = 1 and ε = 10 -5 ). On the top left, we have plotted the distribution function as a function of v and at time t = 1 for x = 25 ∆x. Figure 5 shows that the distribution is not Maxwellian, the collisions are not sufficiently numerous to ensure the relaxation towards the Maxwellian equilibrium. But, on Figure 6 the small Knudsen number yields the convergence towards this Maxwellian equilibrium. Since the Debye length is not small, the quasi-neutral regime is not reached and the density is far from the equilibrium ρ = 1 (see Figures 5 and6, pictures on the top right). The electric fields given by the Poisson equation (28b) and by the reformulated Poisson equation (25b) are identical (see Figures 5 and6, pictures on the bottom left). On the bottom right picture of Figures 5 and6, we can see the time steps, calculated with the CFL condition [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], during all the simulation. Remark that time steps are the same for both schemes and that the time step does not need to be smaller than the Knudsen number which means that both schemes are AP in the fluid limit.

The results for a small Debye length (λ = 10 -5 ) are given on Figure 7. On the bottom right picture, we have plotted the distribution function (just before the code fails) given by the classical scheme with a time step not resolving the plasma period. The scheme is unstable. The other three pictures give the results of the classical scheme with a time step resolving the time step (∆t = 0.9 λ) and those of the AP scheme but with a time step greater than λ as we can see on the top right picture of Figure 7. The AP scheme is stable even for large time steps only resolving the CFL condition [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF] but not the plasma period λ. Furthermore, on the top left picture, we can see the distribution function as a function of v and at the final time t = 8 × 10 -3 = 8000 plasma periods, for x = 25 ∆x. We can see that the classical scheme gives diffusive results due to the fact that it must use small time steps. The diffusion in the velocity space leads to a bad description of the plasma temperature and then to non consistent results. Furthermore, due to this diffusion in the velocity space, the unphysical velocity boundaries ±v max are always attained with the classical discretization and then the classical scheme becomes unstable. This means that the classical number. This means that we only want to observe the macroscopic behaviors of a plasma independently from the small scale dynamics. This is precisely the scope of the present paper and this can be made more clear observing the behaviors of the scheme when the same test is run for different values of λ and for a fixed α here α = 10 -4 . The final time is T = 10. The classical scheme is still stable for λ = 10 -3 but no more for the smaller values even if ∆t resolves the plasma period λ. This is due to the large numerical viscosity of the scheme (see previous test case). Then, we have only plotted the results for our AP scheme. The time step is the same for all values of λ and given by the CFL condition [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF], its value is almost ∆t = 2 × 10 -3 . can observe that the scheme is stable uniformly with respect to λ. Furthermore, the values of λ are small and the scheme is consistent with the quasi-neutral limit since the density is projected on the quasi-neutral state ρ = 1.

Conclusions

In this paper we have discussed a new class of numerical methods which are able to overcome the strong restrictions due to the fast scale dynamics related to collisions and quasi-neutrality in plasmas. In particular we focused on the quasineutrality problem, the collisional dynamics being considered rather simple to treat. The new scheme proposed is able to exactly preserve the quasi-neutrality and the Maxwellian equilibrium constraints for all times and it projects the solution towards the corresponding equilibrium states when the parameter which describes the small scale becomes zero. In addition, we proved the uniform linear stability for all values of the Debye length and the projection to the equilibrium state after one time step when this length is set to zero. We tested the performances and the stability of the scheme with several numerical simulations using different values of the scaling parameters. The method presented permits to avoid the resolution of the fast microscopic dynamics and allows to directly jump to the analysis of the macroscopic phenomena without forcing the numerical parameters like the time and the space steps to be constrained by these small scales dynamics.

We are actually working on the extension of the present technique to high order in time and space schemes which share the same properties of the one presented in this paper. This will permit, for instance, to be more accurate in studying the phenomena of relaxation towards the steady states macroscopic regimes. We also would like to study more physical multi-dimensional phenomena, to treat more realistic collisions models, such as Landau-Fokker-Planck operators and to add magnetic fields to the problem. ). The slope of the line is -0.85 (given by the Landau damping rate calculated in [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF]). Right: density as a function of x and at time t = 5. ). The slope of the line is -0.1 (given by the Landau damping rate calculated in [START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF]). Right: density as a function of x and at time t = 10. The time step, ∆t ≈ 4 × 10 -3 , is given by the C.F.L. constraint (29).

Figure 2 :

 2 Figure 2: Numerical value of y 0 solution of equation (27)

Figure 5 :

 5 Figure 5: Periodic perturbation test case for ∆t < ε = λ = 1, comparison of the classical and AP schemes. Top left: distribution function as a function of v for x = 25 ∆x and t = 1. Top right: time steps as functions of time. Bottom: density and electric field as functions of x at t = 1.

Figure 6 :

 6 Figure 6: Periodic perturbation test case for ε = 10 -5 < ∆t < λ = 1, comparison of the classical and AP schemes. Top left: distribution function as a function of v for x = 25 ∆x and t = 1. Top right: time steps as functions of time. Bottom: density and electric field as functions of x at t = 1.

Figure 7 :

 7 Figure 7: Periodic perturbation test case for ε = λ = 10 -5 , comparison of the classical and AP schemes. Bottom right picture: unstable distribution function given by the classical scheme when ∆t > λ. Left pictures: distribution function as a function of v for x = 25 ∆x and density as a function of x both at time t = 1 and for time steps given on the top right picture.

Figure 8 :

 8 Figure8: Landau damping test case: λ = 1, α = 10 -2 . Left: time evolution of log( E 2 ). The slope of the line is -0.85 (given by the Landau damping rate calculated in[START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF]). Right: density as a function of x and at time t = 5.

Figure 9 :

 9 Figure9: Landau damping test case: λ = 0.1, α = 10 -2 . Left: time evolution of log( E 2 ). The slope of the line is -0.1 (given by the Landau damping rate calculated in[START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF]). Right: density as a function of x and at time t = 10. The time step, ∆t ≈ 4 × 10 -3 , is given by the C.F.L. constraint[START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF].

Figure 10 :Figure 11 :

 1011 Figure10: Landau damping test case: λ = 0.1, α = 10 -2 . Left: time evolution of log( E 2 ). The slope of the line is -0.1 (given by the Landau damping rate calculated in[START_REF] Belaouar | An asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF]). Right: density as a function of x and at time t = 10. The time step is constant: ∆t = 10 -4 .

the authors were supported by the French ANR project Moonrise. N. Crouseilles and M.H. Vignal were supported by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03. N. Crouseilles was supported by the ERC starting grant project GEOPARDI.

scheme can not be stabilized for large times whereas the AP scheme can.

Linear Landau damping

In this Section, we only consider the Vlasov-Poisson system and we neglect the collisions, then we consider system (4) with ν = 0 which corresponds to the Vlasov-Poisson system.

We study the one dimensional linear Landau damping test case. This test case consists in a perturbation of a stationary quasi-neutral plasma described by the Vlasov-Poisson equation. We initialize the problem with

where 0 < α 1 will be specified for each simulations. The space domain is x ∈ [0, 2 π]. The problem is supplemented by periodic boundary conditions in space for the distribution function f and by homogeneous Dirichlet boundary conditions for the electric potential, homogeneous Dirichlet boundary conditions are used as well for the distribution f in the velocity direction. The velocity space is truncated at v max = ±12 and discretized by N v = 128 points while the spatial direction is discretized by N x = 128 points. The time step is limited by the CFL condition [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF].

For all simulations, we plot the time evolution of log( E 2 ), and the density ρ as a function of x for a given time. We recall that the potential energy is given by λ/2 E 2 (x) dx = λ/2 E 2 2 . First, we consider λ = 1 and α = 10 -2 while the final time is T = 10. On Figure 8, we can see that the AP scheme captures the solution in term of the Landau damping as well as the classical explicit scheme does. On the right picture of Figure 8, we observe that the density ρ(t = 5, x) is far from the quasi-neutral equilibrium for both schemes.

In the second test, we consider λ = 0.1, α = 10 -2 and a final time T = 10. On Figure 9, we observe that the AP scheme presents a stronger damping rate and a decay to the equilibrium for density compared to which is faster than the one computed with the classical scheme. This is natural and expected because the AP scheme is constructed to be consistent with the quasi-neutral limit and thus its damping behaviors are larger than those of the classical scheme which can not capture the asymptotic quasi-neutral regime. Of course if one is interested in the microscopic dynamics, which is described by the small scale, one must be in a resolved regime, this means that one must take smaller time steps. On Figure 10, we plot the results for a time step given by 10 -4 (almost 40 times smaller than that used for the results of Figure 9). We observe that the AP scheme converges towards the exact solution.

On Figure 11, we consider different small values of λ for α = 10 -4 . We recall that we are interested in constructing a numerical scheme which projects the solution to the quasi-neutral state with a stability condition which does not depend on the small scales induced by the Debye length and the Knudsen