KINETIC/FLUID MICRO-MACRO NUMERICAL SCHEMES FOR VLASOV-POISSON-BGK EQUATION USING PARTICLES

. This work is devoted to the numerical simulation of the Vlasov equation in the ﬂuid limit using particles. To that purpose, we ﬁrst perform a micro-macro decomposition as in [3] where asymptotic preserving schemes have been derived in the ﬂuid limit. In [3], a uniform grid was used to approximate both the micro and the macro part of the full distribution function. Here, we modify this approach by using a particle approximation for the kinetic (micro) part, the ﬂuid (macro) part being always discretized by standard ﬁnite volume schemes. There are many advantages in doing so: ( i ) the so-obtained scheme presents a much less level of noise compared to the standard particle method; ( ii ) the computational cost of the micro-macro model is reduced in the ﬂuid regime since a small number of particles is needed for the micro part; ( iii ) the scheme is asymptotic preserving in the sense that it is consistent with the kinetic equation in the rareﬁed regime and it degenerates into a uniformly (with respect to the Knudsen number) consistent (and deterministic) approximation of the limiting equation in the ﬂuid regime.


Mohammed Lemou
Abstract. This work is devoted to the numerical simulation of the Vlasov equation in the fluid limit using particles. To that purpose, we first perform a micro-macro decomposition as in [3] where asymptotic preserving schemes have been derived in the fluid limit. In [3], a uniform grid was used to approximate both the micro and the macro part of the full distribution function. Here, we modify this approach by using a particle approximation for the kinetic (micro) part, the fluid (macro) part being always discretized by standard finite volume schemes. There are many advantages in doing so: (i) the so-obtained scheme presents a much less level of noise compared to the standard particle method; (ii) the computational cost of the micro-macro model is reduced in the fluid regime since a small number of particles is needed for the micro part; (iii) the scheme is asymptotic preserving in the sense that it is consistent with the kinetic equation in the rarefied regime and it degenerates into a uniformly (with respect to the Knudsen number) consistent (and deterministic) approximation of the limiting equation in the fluid regime.

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
use the so-called mean free path, which is the relative distance (or time of flight) made by a particle between two consecutive collisions. Due to their high dimensionality and their multiscale character, the numerical resolution of kinetic models for plasmas is challenging. For instance, due to their prohibitive numerical complexity, phase space grid methods cannot be used to accurately describe the full 6 dimensional case. On the other side, particles methods are extensively used in real situations due to their low computational cost [2]. However, these particle methods are affected by an important noise due to their probabilistic character, and standard time discretizations can be constrained by a severe stability condition depending on a small parameter such as the mean free path.
The main goal of this work is to design a particle-based Asymptotic Preserving (AP) scheme for a Vlasov-Poisson-BGK equation in the hydrodynamic limit. This concept of AP scheme has been introduced in [21]. Using the micro-macro approach (see [23,3,9]), we decompose the distribution function into an equilibrium part and a fluctuating part. The original kinetic equation is then equivalently reformulated into a coupled system where the macro part is an Euler type model and the micro part is a Vlasov type equation. Our general strategy is to discretize the macro part (of Euler type) on a spatial grid whereas the kinetic part is approximated using a particle method. We faced to several difficulties due to the particle nature of our approach. First, the macro part of the coupled system involves fluxes of the kinetic part and therefore the particle approximation of this kinetic part should be projected on a spatial grid to be used in the deterministic computation of the macro part. Second, the property of the kinetic part to have vanishing first moments should be ensured along the simulation to be coherent with the micro-macro structure of the coupled system. Indeed, the direct computation of the moments of the kinetic part does not ensure this property and this violates the important conservation properties of the system, and in practice is a source of numerical fluctuation and undesirable noises. Therefore, an adjustment procedure is needed. Moreover, the particle resolution of the kinetic part needs a splitting method between the transport part and the various source terms. Finally, a suitable semi-implicit strategy enables to overcome the difficulty induced by the stiff source term and leads to desired Asymptotic Preserving property.
Consequently, this strategy is Asymptotic Preserving and allows the use of a time step ∆t which is independent of the mean path (see also [6,18]). Note that this time step is also independent of the usual CFL transport condition of the kinetic part since we use here a particle approach on this part; however the time step remains constrained by the hyperbolic structure of the macro part.
In addition, the so-obtained numerical scheme satisfies the following interesting properties: (i) the so-obtained scheme presents a much less level of noise compared to the standard particle method; (ii) the computational cost of the micro-macro model is reduced in the fluid regime since a small number of particles is needed for the micro part; (iii) the scheme is asymptotic preserving in the sense that it is consistent with the kinetic equation in the rarefied regime and it degenerates into a uniformly (with respect to the Knudsen number) consistent (and deterministic) approximation of the limiting equation in the fluid regime.
There is an important literature dealing with the construction of suitable AP schemes for kinetic equations in various contexts. We mention for instance works based on domain decompositions, separating the macroscopic (fluid) domain from the microscopic (kinetic) one (see [8,14]). This kind of methods faces a natural KINETIC/FLUID SCHEMES FOR VLASOV USING PARTICLES 3 difficulty which is linked to the handling of the interface between the two domains. There is another kind of AP schemes for kinetic equations, which are based on the use of time relaxed techniques where the Boltzmann collision operator is discretized by a spectral or a Monte-Carlo method (see [28,27,19]). Other techniques have also been developed to design multiscale numerical methods which are based on splitting strategy [6], penalization procedure [17,18] or micro-macro decomposition [23,3,9]. We choose to follow the micro-macro strategy since it is a systematic method in the sense that it can be applied to different asymptotics (diffusion, fluid, high-field, ...).
Most of these approaches are performed on a phase space mesh and the numerical computation has a constant cost with respect to the Knudsen number. Despite of this AP property, the computational cost of this method still needs to be reduced in the fluid regime, since in this regime, one does not need a refined grid in the velocity direction and only few velocity points may be sufficient. To remedy this problem, our strategy here is to couple a particle method with an AP strategy based on a micro-macro decomposition. Indeed, our strategy uses a particle method for the micro (kinetic) part of the decomposition and this allows the use of few particles in the fluid regime. To do that, we adopt an approach similar to that proposed in [11,16,26,20,15] in order to obtain a particle discretization of the micro unknown in the coupled model. An adapted semi-implicit discretization enables to design a particle AP scheme in the fluid limit. Our approach finally bears some similarities with the moment guided particles method [12,13] or with the delta-f method [4,5]. But here, the numerical scheme enjoys the AP property and the particle approximation is only used on the micro part which allows to reduce its cost in the fluid regime. Additionally, a coupling with the Poisson equation is also considered in our model.
To wit, we consider a distribution function f (t, x, v) ≥ 0 of electrons in phase space, and E(t, x) the self consistent electric field. The fluid scaling for the Vlasov-BGK equation reads ∂f ∂t coupled with the Poisson equation where x ∈ [0, L] and v ∈ R are the phase space independent variables and t ≥ 0 the time. The collision operator Q(f ) is the BGK collision operator where M (U ) is the Maxwellian associated with f : and U is the vector of the first moments of f (density, momentum and energy)

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
The functions f and E are submitted to the following periodic condition In order to get a well-posed problem, a zero-mean electrostatic condition has to be added, together with an initial condition The parameter ε denotes the Knudsen number. When ε is small, the mean free path becomes small compared to the size of the domain, the kinetic model degenerates (at least formally) into Euler-Poisson equations or Navier-Stokes-Poisson equations, which are the relevant models in this regime. However, when ε is of order one, these models can not correctly describe the plasma.
As in [3], we deal with the simple BGK model for which the micro-macro model is used to design an AP scheme in the fluid limit. As said before, a self-consistent electric field is considered here. An important input in this work is the fact that the kinetic part of the micro-macro model is discretized by using particle method, whereas a grid of the phase space is used in [3]. In this way, we construct a scheme with the following two properties: (i) the usual noise which is observed in standard PIC method is strongly reduced by the micro-macro decomposition, (ii) the PIC method is used only for the micro part and the number of particles can be very small in the fluid regime. This avoids unnecessary calculations in this regime and reduces the computational cost. We emphasize that despite the noisy character and the slow convergence of the standard PIC method, they are extensively used in real situations (requiring high dimensions) due to their large flexibility and low computational cost. We refer to [2] for physical applications and to [7] for a mathematical analysis. In our context, one of the main difficulties is to maintain the micro-macro structure along the time evolution; indeed, the first three moments in velocity of the kinetic part must be zero for all time and a suitable numerical scheme is constructed here to guarantee this property. In the spirit of the matching procedure proposed in [12], we proposed an additional step of the algorithm which ensures that the moments of the kinetic part of the decomposition are strictly zero.
The remainder of the paper is organized as follows. Some basic properties of the Vlasov-Poisson-BGK model and its fluid limit are recalled. The derivation of the micro-macro model is then explained. Details of the numerical method are given in section 4. Finally, numerical results are proposed to illustrate the efficiency of the method. 2. BGK model and its fluid approximation. Now we briefly describe the wellknown conservation laws of (1) and its asymptotic models when ε goes to zero. Multiplying (1) by m(v) and integrating with respect to v yields This is equivalent to the following system

hal-00728875, version 1 -6 Sep 2012
When ε goes to zero in (1), the distribution function tends to a local Maxwellian M (U ) given by (4). The previous system can then be closed, and the heat flux v 3 f = v 3 M (U ) can be expressed as a function of U Using a Chapman-Enskog expansion, corrective terms of order ε can be derived to this model which leads to the usual compressible Navier-Stokes equations for plasmas. More precisely, in our one dimensional context, the two first equations (density and momentum) are unchanged whereas the energy equation involves the following source term −ε∂ x (κ∂ x T ), where the heat conductivity κ = (3/2)ρT depends on U (see [1]).
3. Derivation of the micro-macro model and the asymptotic limit. This section is devoted to the derivation of the micro-macro model starting from the Vlasov-BGK equation following the lines in [3]. The only difference with the decomposition in [3] is the presence of a self-consistent electric field which we have to incorporate in the decomposition.
3.1. Derivation of the micro-macro model. We first write f according to the following decomposition (when no confusing is possible, we will use the notation with m(v) = (1, v, v 2 /2) T and where the macroscopic variables U (t, x) = (ρ(t, x), (ρu)(t, x), (ρu 2 /2 + (1/2)ρT )(t, x)) are the first three moments of f In particular mgdv = 0. Let T be the transport operator T f = v∂ x f + E∂ v f , then the Vlasov equation (1) writes We denote by Π M the orthogonal projection in L 2 (M −1 dv) endowed with the weighted scalar product (ϕ, ψ) M = ϕψM −1 onto the following space with N (L Q ) the null space of the linearized operator L Q of Q. The explicit expression of this projection operator is given by (see [3]) (9) This projection will be used to derive from (8) a macro equation on M (or equivalently on U ) and a micro equation on g. Applying (I − Π M ) to (8) gives

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
Since (I − Π M )(∂ t M ) = Π M (g) = Π M (∂ t g) = 0 (Lemma 3.1 of [3]), we get the micro part Applying now Π M to (8) leads to Since T M = v∂ x M + E∂ v M , we obtain the macro part of the micro-macro system, where the term F (U ) (which corresponds to the usual Euler fluxes) and the source term S(U ) are given by Finally, the micro-macro model of unknown (g, U, E) can be written as follows where Π M is defined in (9), U = (ρ, ρu, ρu 2 /2+ρT /2) T g = v∂ x g +E∂ v g, the Euler fluxes F (U ) and the source term S(U ) are given by (10), · denotes the integration in v and m(v) = (1, v, |v| 2 /2) T . It is formally clear that system (11)-(12)-(13) is equivalent to the original kinetic equation. This statement is summarized in the following proposition.

3.2.
Chapman-Enskog expansion. We briefly recall how the limiting model (ε → 0) is formally obtained from (1). Since the micro-macro model is equivalent to the original kinetic model, the Chapman-Enskog procedure will be applied to (11)-(12)- (13). From (11), we clearly have g = O(ε) when ε → 0, and then we have which, injected in (12), gives The last term corresponds to the Navier-Stokes correction terms (see [3]) and can be computed to get 4. Numerical approximation. In this section, we introduce our numerical method which is based on a suitable discretization of the micro-macro model (11)-(12)- (13).
The main difference of our approach compared to [3] is the fact that the micro part is solved by a particle method whereas a phase space grid is used in [3]. Obviously, such a method is intended to be faster than an Eulerian method (especially in high dimensions). Moreover, when one deals with hydrodynamic regimes, the function g is small (of order ε) and it becomes unnecessary to keep a refined grid of the phase space. In this case, a coarse discretization of the micro part (with few points) may be sufficient, and a particle method is most suited to this context. As said before, one difficulty of our approach is to maintain the structure of the micro-macro decomposition along the time: mg(t) = 0, ∀t ≥ 0, at the discrete level. This property is of course satisfied by the conitnuous micro-macro model. Indeed, it is observed numerically that a violation of this property generates too much noises in the numerical simulations. This noise is instead strongly reduced when the structure is exactly preserved during the time evolution. This step shares some similarities with the matching procedure of [12] or the delta-f method [4,5]. The approach here is made systematic by the use of a numerical projection which is directly derived from the operator projection Π M defined above. Roughly speaking, one iteration of the present method can be summarized as follows: • solve the micro part (11) using a particle method • suitably modify the obtained function g to ensure the zero-moments property mg = 0 • solve the macro part (12) with a finite volume method where particles are used to evaluate vmg . The particle method is attractive in the fluid regime since few particles will be sufficient to approximate g. Furthermore, in kinetic regimes, we observe that, for a given number of particles, the present method enables to reduce statistical fluctuations which may be observed in many particle approaches, as the so-called δf method (see [4,5]).
In the sequel, we will consider a uniform spatial grid x i = i∆x, i = 0, . . . , N x , with ∆x = L x /N x , N x is the number of points and L x the size of the spatial domain. Periodic boundary conditions are imposed. The time step is ∆t (t n = n∆t, n ≥ 0), the velocity direction is approximated by an interval of length L v and N part denotes the number of macro-particles.
Hereafter, the main steps of the algorithm are presented.

4.
1.1. Generalities on particle methods. The classical approach to solve collisional kinetic equation using particles method is the time splitting between the transport part ∂ t g + T g = 0,

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
with T g = v∂ x g + E∂ v g and the collision (or source) part where S(g) denotes a general source term (such as a collision operator for example) which depends on t, x and v.
In particle methods, the distribution function g is approximated by a finite set of particles where x k (t) represents the position, v k (t) the velocity and ω k (t) (k = 1, . . . , N part ) the weight of the k-th particle (see [29,20,26,16,15,11]). In particular, ω k (t) and the function g are related through During the first step of the splitting (transport step), the weights ω k are constant in time and the characteristics have to be solveḋ This step corresponds to the standard Particle-In-Cell (PIC) method (see [2]). We move particles thanks to the equations of motion, whereas the electric field E is obtained by solving the Poisson equation. A standard pusher (Euler for example) can be used: knowing x n k , v n k (which are some approximations of x k (t n ) and v k (t n )) and E (t n , x n k ), we compute x n+1 Obviously, more accurate pushers can be used (such as the second order Verlet scheme, see [30]). In general, E(t n , x) is known on the uniform spatial grid (x i ) i (which makes easier its resolution through the Poisson equation using FFT) so that an interpolation is necessary to compute E(t n , x n k ). When one deals with a source term S, the second step of the splitting enables to modify the weights ω k (t) through the following ordinary differential equatioṅ where s k (t) is the weight associated to the source term S. Here we use the following definition that links an approximation of the function at (x k , v k ) with its corresponding weight Any solver can be applied to numerically resolve (19). A simple example is the forward Euler scheme ω n+1 k = ω n k + ∆t s n k . Note that L v may depend of time since v max can evolve from t n to t n+1 .
We also detail the computation of the moments of the distribution function g given by (16), on the uniform grid (x i ) i . To do that, the particles are regularized using a convolution function as the B-splines for example (see [2]). Starting from the particles representation of g (16) the j-th moment in v of g is computed at position where x k , v k , ω k denote the position, the velocity and the weight of the k-th particle, (x i ) i the uniform spatial grid and B ℓ ≥ 0 is a B-spline function of order ℓ: Numerical tests have been performed up to first order (ℓ = 1). Now, we describe our strategy to solve the micro-macro system (1)-(2)-(3). As said before, this is done by following three steps: numerical resolution of (11) using particles, matching the moments of g to zero and numerical resolution of (12) using a finite volume method.

4.1.2.
Numerical resolution of (11) using particles. In this step, we propose a splitting procedure for the numerical resolution of (11) by means of a particle method.
Following the previous part, we consider the following splitting The first part of this splitting is performed following a numerical solver for (16) such as (18). The second part needs more attention. It can be reduced to the resolution of the following ordinary differential equatioṅ where α k (t) is the weight associated to (I − Π M )(T M ) We refer to the next subsection for the details of the numerical computation of Π M .
To design an Asymptotic Preserving solver, it is necessary to consider an implicit discretization of the stiff term ω k (t)/ε in (22). Then, we consider the following discretization of (22) which ensures the stability with respect to ε for any fixed ∆t. Now we say few words on the computation of α k (t) and β k (t). As said before, we consider the following relation (24) We see in Appendix A the expression of (I − Π M )T M and Π M (T g) and observe that they depend on x and v and have the following polynomial form in v: where the coefficients a ℓ = a ℓ U, ∂ x U, ∂ x v 3 g depend on the macro unknown U , its spatial derivative ∂ x U and on ∂ x v 3 g (see (32) and (29) in Appendix A). The Maxwellian M only depends on U and its form is explicit in velocity. Consequently, once (I −Π M )T M and Π M (T g) are computed, the weights α k and β k can be obtained using (24). Considering that U is known at the spatial grid points (x i ) i and that (x k , v k , ω k ) is known, we obtain α k and β k as follows 10

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
• compute the spatial derivatives of U on spatial grid points (x i ) i using finite differences formula, • compute v 3 g on spatial grid points (x i ) i using (20), and compute its spatial derivative using finite differences formula, • following (24), evaluate the so-obtained functions (polynomial in v times Maxwellian) at (x k , v k ); a linear interpolation is used to get U (x k ) using U on the spatial grid points (x i ) i .

4.2.
Projection step: matching the moments of g to zero. In this subsection, we present the discrete version of Π M which will ensure in particular that, at each time step, mg = 0. In fact, nothing guarantees that this property is satisfied at the discrete level, since mT g = 0 and the set of functions g such that mg = 0 is not stable by the action of the operator T .
To compute the projection Π M of a g n+1 , we seek a function in the kernel of the linearized collision operator N (L Q ) which has the same first three moments of g n+1 . Such a function, which we call h(x, v), has the following form and satisfies mh | x i = mg n+1 | x i . The main idea is then to determine λ(x) ∈ R 3 such that mh | x i = mg n+1 | x i at the discrete level (three unknown for three constraints). Hence, the micro unknown is replaced by (I − Π M )g n+1 = g n+1 − h which satisfies by construction m(g n+1 − h) | x i = 0, ∀i = 1, . . . , N x .
In practice, the procedure is applied to the weights. Denoting by γ k the weight associated to the function h: the weights ω n+1 k of g n+1 are then replaced by ω n+1 k − γ k so that, at the discrete level we ensure that the first three moments of these new weights are equal to zero, as desired. It corresponds somehow to the application of a discrete approximation of the operator (I − Π M ) which is of course consistent with the continuous model.
We now detail the computation of the function h (which is the discrete version of Π M (g n+1 )). We will denote by p k the weight of the Maxwellian M , according to the relation (24): First, we expand λ(x) on the basis of B-splines of degree ℓ given by (21):

hal-00728875, version 1 -6 Sep 2012
Then, the moments of h at x i are Then, we have to look for λ ∈ R 3Nx solution to the linear system U g = Aλ where λ ∈ R 3Nx is the vector whose components are λ j ∈ R 3 , j = 1, . . . , N x , and A is a so that we have to solve ∀i = 1, . . . , N x the linear system U g (x i ) = Nx j=1 A i,j λ j . We detail in the sequel the first orders which have been tested ℓ = 0, 1. We recall that with B 0 given by (21), so that the moments of h at x i are Then, the equality U g (x i ) = mh | xi is equivalent to the following 3 × 3 linear system In this case, the (3N x × 3N x ) system is decoupled so that we have to deal only with N x (3 × 3) linear systems. Once the coefficients λ i ∈ R 3 , ∀i = 1, . . . , N x are determined, the weights γ k of h are computed as usual using where j k is such that the k-th particle x k satisfies |x k − x j k | < ∆x/2. The weights ω n+1 k of g n+1 are then replaced by ω n+1 k − γ k . All these computations are done at the discrete level which allow to get exact discrete conservations up the machine precision.

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
Order 1. We want to extend the procedure to the first order B-spline function B 1 We now compute the moments of h at x i as and if x k ∈ I i−1 , we have Hence, we finally have for the moments of h, where A i and B i are given by (for i = 1, . . . , N x )

hal-00728875, version 1 -6 Sep 2012
and M j,i , N j,i , j = 0, . . . , 4, i = 1, . . . , N x , We can write this relation as a (3N x × 3N where 0 3 is the zero 3 × 3 matrix and where periodic conditions have been used. Finally, the weights ω n+1 k of g n+1 are replaced by ω n+1 Remark 1. This strategy can be generalized quite easily to higher order B-spline functions. The more the order is high, the more the matrix A involved in the linear system will be large and full which can make its inversion expensive. However, the matrix remains symmetric and has a block-diagonal structure which can be optimized by using a specific linear solver so that this step does not cost too much compared to the rest of the algorithm. (12). For the macro part, we use standard finite volume methods on the uniform mesh (x i ) i , such as Lax-Friedrich or Rusanov scheme. Obviously, many efficient numerical methods from the literature can be applied (see [24]). We have to solve the macro equation (12):

Numerical resolution of
where F (U ) are the Euler fluxes and S(U ) = (0, ρE, ρuE) which can be written where F n i+1/2 is chosen to be the Rusanov flux, that is: where a i+1/2 is the maximum of the absolute value of the eigenvalues of the Jacobian of F over cell i and (i + 1). The source termS is approximated bỹ where x i+1/2 denotes the middle of the cell [x i , x i+1 ], ∀i = 0, . . . , N x − 1. The terms which involves g are nothing else but moments of g whose computation is performed as in (20).

4.4.
Asymptotic preserving property. In this subsection, we prove (formally) that the numerical scheme presented above enjoys the asymptotic preserving property. Indeed, solving the micro part (23) leads to where α n k denotes the weight of (I − Π M n )(T M n ) and β n k the weight of Π M (T g n ). When ε is small, we directly have from (26) that ω n+1 k = O(ε) and the numerical scheme degenerates into a consistent discretization of the Euler equations. Moreover, if we consider the first order in ε for ω n+1 k , we have from (26) Since the position and velocity of the particles do not move during this step, we obtain after this step which, in terms of the distribution function, means that This is exactly the non equilibrium part obtained in section 3.2. Then, the moments of vg n+1 have to be computed in the macro equation (25). This is done following (20). We then obtaiñ The last term is a consistent approximation of the corrective terms in the Navier-Stokes equation which ensures that the present scheme enjoys the desired asymptotic preserving property for the Navier-Stokes equation.

Remark 2.
As suggested in [22,9], it is possible to extract the diffusion term coming from (I − Π M )(T M ) in order to consider it in the macro part. Injecting the expression (23) of g n+1 in the kinetic fluxes of the macro equation (25), we get where the term N i+1/2 is the Navier-Stokes diffusion The consequence of this operation is that the Navier-Stokes diffusion term can be implicited as in [22,9] so that the so-obtained scheme enjoys the AP properties and the associated asymptotic numerical scheme involves implicit diffusion Navier-Stokes terms.

4.5.
Algorithm. The global algorithm is then the following: 1. push forward the particles using the scheme (18), 2. compute, on the spatial grid, using deposition formula (20) the quantities T M , Π M (T M ) and Π M (T g), 3. compute the new weights using (23), 4. replace g by (I − Π M )(g) to ensure that the moments of g are zero (see subsection 4.2), 5. compute v 3 g using (20), 6. advance the macro equation using (25). The numerical cost of this algorithm is of the same order as a particle solver for the Vlasov-Poisson-BGK model (1)-(2)-(3). Indeed, the additional steps involved by the micro-macro model have a cost of the order of the number of particles N part . In practice, the numerical cost of the micro-macro model is about twice the cost of the standard PIC method on the Vlasov-BGK model. However, our strategy to use a particle method for the kinetic part (equation on g) allows to capture the fluid regime with few particles, as we will see in the next section.

5.
Numerical results. We present here some numerical results obtained with the micro-macro model (11)-(12)-(13) using the algorithm presented above. We are interested in the simulation of classical plasma test cases. When ε is of order 1, the micro-macro model (called MiMa) will be compared to the classical PIC discretization of the Vlasov-BGK model (denoted by PIC-BGK), whereas for small values of ε, it will be compared to the Navier-Stokes or Euler models (referred to as NS and Euler). These two last models are approximated using a finite volume method with a Rusanov approximation of the fluxes as presented in subsection 4.3. We also had some results obtained by a deterministic solver of the Vlasov-Poisson-BGK equation. 5.1. Test 1: linear Landau damping. We first consider a Landau damping test case in which the amplitude of the perturbation is small. In this case, the nonlinear system is close to its linearized form. The initial condition is then (with α = 0.01) with the wave number k = 0.5. For the micro-macro model, the initial condition is U (t = 0, x) = (1 + α cos(kx), 0, 1 + α cos(kx)) and g(t = 0, x, v) = 0.
The numerical parameters are considered as follows: • Euler: N x = 128, ∆t = C∆x, with C = 0.4. Indeed, the PIC method does not induce any restriction on the time step but an explicit treatment of the weight equation (collisional part) induces a condition ∆t = Cε for the PIC-BGK model. On the other side, the MiMa model is only restricted on the macroscopic stability condition due to the Euler solver: the maximum velocity being |u| + √ 3T , the CFL condition then writes ∆t < ∆x/(|u| + √ 3T ). We take several values for the number of particles in the PIC-BGK and MiMa methods to check its influence in the different regimes: fluid, intermediate and non-collisional regimes.
For this test, we are interested in the time evolution of the electric energy whose oscillations are known to damp exponentially in time, according to the collisionless Landau theory. We represent the evolution of the electric energy (27) in a semi-logarithmic scale. Following the collisionless Landau theory, this quantity presents oscillations whose amplitude decreases with a linear rate. When a collision operator is considered, this remains true, even if the rate now depends on ε. More precisely, as ε goes to zero, this rate is expected to converge to zero since the fluid regime does not capture the Landau damping. Figures 1 and 2 illustrate this fact by comparing MiMa (micromacro model) for different values of ε to the Euler equation. The AP property is also confirmed in this figure since with a fixed ∆t, the results of the micro-macro models become closer to the Euler model as ε goes to zero. In particular, when ε = 10 −4 (Figure 2), the fluid regime is well simulated since the two curves are nearly the same.
In Figure 3, we compare the results of MiMa and PIC-BGK against NS (Navier-Stokes equations) for ε = 0.1. The number of particles has been chosen so that the numerical convergence is ensured for MiMa. We observe that the three models produce the same results up to t ≈ 10; beyond this time, PIC-BGK does not provide the correct behavior because of its inherent numerical noise. On the contrary, MiMa and NS are very close up to the end of the simulation. In the intermediate regime (ε = 0.1), MiMa has also a very good behavior compared to PIC-BGK.
We are now interested in larger values of ε. We will compare the micro-macro model to a PIC discretization of the Vlasov-BGK model. First, on Figures 4 we compare the electric energy obtained with MiMa and PIC-BGK for ε = 1 and ε = 10. We observe that PIC-BGK cannot be accurate enough when it deals with values which are close to the level of numerical noise (after t ≈ 15). The results of MiMa are really good in this case since it is able to reproduce the correct behaviour of the electric energy (the slope is equal to −0.16 which is very close to the theoritical value −0.1533 see [30]) for very long time and up to round-off errors, using the same number of particles as PIC-BGK. Now we want to know if it is possible to decrease the level of numerical noise of PIC-BGK simulations by adding particles. On Figure 5, we show that it is true, but PIC-BGK is not as accurate as MiMa, even if more particles are considered. One explanation is the following: MiMa uses particles to describe only the microscopic part g and not all the distribution function f ; then, the numerical noise only affects g (which remains quite small in this test) for MiMa whereas all the distribution function f is affected by the noise. It is difficult to quantify the improvement of adding particles, but looking at Figure 5, we see that, even with a large number of particles, the PIC-BGK discretization is far from reaching the same low level of noise as MiMa.
We are now interested in spatial dependent diagnostics. In Figure 6, we plot the charge density (1 − ρ) at two fixed times (t = 0.2 and t = 0.4), as a function of x, for MiMa and PIC-BGK with ε = 1. As already observed before, when the number of particles is fixed, the numerical noise is strongly reduced using the micro-macro decomposition compared to the standard PIC-BGK approach.
On Figure 7, we plot the full distribution function f = M + g. We compare the one given by MiMa to the one obtained with the PIC-BGK approach, at time t = 20 and for ε = 1. Here, it is an illustration of the numerical noise of the PIC-BGK approach compared to the MiMa one. The numerical noise arising in the MiMa approach is very small since it only concerns the fluctuation g. This is emphasized by Figure 8 in which we plot the difference of the distribution functions obtained by MiMa and PIC-BGK.
We have also looked at conserved quantities. We have first verified that the total mass f dvdx and the total momentum vf dvdx are conserved exactly. The total energy [ v 2 f dvdx+ E 2 dx] is not exactly preserved but variations of about 0.1% are observed for MiMa, for all ε > 0.
Finally, we look at the behavior of the heat flux as a function of ε. Indeed, we have seen in subsection 4.4 that the semi-implicit discretization of the micro-macro model allows to recover the Navier-Stokes asymptotics. To check this point, we compare the heat flux given by the Navier-Stokes equation (−3/2)ρT ∂ x T to the heat flux given by the micro-macro model (1/ε) (v − u) 3 f dv. This last quantity can be simplified into (since (v − u) 3 M dv = 0) (1/ε) (v − u) 3 gdv; moreover, since mg = 0, we finally have to compute the third order moment of g: (1/ε) v 3 gdv using (20). The difference of these two quantities has to be of the order of ε. In Figure 9, the difference (in L 2 norm) of the two heat fluxes (obtained with MiMa and NS) is plotted as a function of ε, at t = 1; different numbers of particles are used to observe the influence of N part . The good order is recovered since a linear behavior is observed. This is illustrated by Figures 10 where the two heat fluxes (for NS and MiMa) are plotted at time t = 1 as functions of x. We observe that when ε decreases, the two heat fluxes become closer.

Test 2.
This subsection is devoted to a second test case, in which the initial distribution function does not belong to the kernel of the linearized collision operator. We take α = 0.05 in with the same wave number k = 0.5. For the micro-macro model, we have the following initial condition for U (macro part) U (t = 0, x) = (1 + α cos(kx), 0, 5 (1 + α cos(kx)) , and for the micro part g (see computations in Appendix B) For the Euler and Navier-Stokes equations, we consider U (t = 0, x) = (1 + α cos(kx), 0, 5 (1 + α cos(kx))).
We are interested in the same diagnostics as in the previous test case. The same numerical parameters are also considered; the Euler equation is now solved using a CFL number C = 0.1. On Figures 11 and 12, we represent the electric energy given by (27) for MiMa and Euler: we study the convergence of the micro-macro model to the Euler equation when ε goes to zero. The AP property is emphasized in Figures 13 in which the L 1 norm of g is plotted as a function of time. For the L 1 norm, we choose g L 1 = ( k |ω k |)/(L x L v ). As in [31], we observe that g becomes of order ε even if it is not at time t = 0 (note that g(t = 0) L 1 ≈ 5 × 10 −2 ). The AP property is satisfied for this initial condition which is not close to equilibrium, even if a quite low number of particles is used.
When larger values of ε are considered, we see on Figure 14 that MiMa is able to give correct results compared to PIC-BGK which suffers from its inherent noise since the electric energy cannot damp towards zero.
As previously, we look at the reconstructed distribution functions for PIC-BGK and MiMa on Figures 15, at time t = 10. Even if g is not small, we can observe that the noise is very smaller than in the PIC-BGK function. This could be an explanation for the quite good behavior of MiMa compared to PIC-BGK. Finally, we verified that the total mass and momentum are exactly preserved. It is not the case for the total energy [ v 2 f dvdx + E 2 dx], but as in the previous case, it is also well preserved (about 0.1% for all ε > 0).

5.3.
Computational cost. One of the main objective of this work is to compute at a low numerical cost the fluid regime with a kinetic model. On the one side, the micro-macro approach enables to design an asymptotic preserving numerical scheme so that the numerical cost is independent of ε. On the other side, the use of particles allows us to decrease the computational cost since the amplitude of g decreases as ε becomes smaller. Indeed, the use of a refined grid is not necessary for small values of ε and our claim is that only few particles are needed to compute the small quantity g.
In this subsection, we study the number of particles when ε varies. For each given ε, we associate the number of particles N part such that the convergence is reached. On Figure 16, we observe that the convergence is obtained for N part = 5 × 10 3 when ε = 10 −2 and N part = 500 is sufficient when ε = 10 −4 .
From our numerical observations, for ε in the range [0.001, 0.1] we may extract the following empirical law N part ≈ C 1 exp(C 2 ε C3 ) for N x = 128. This is shown in Figure 17, with C 1 ≈ 2000, C 2 ≈ 15.5 and C 3 ≈ 0.6.
Since the complexity of MiMa is O(N part ), N part being always greater than N x , the gain in terms of computational cost is greatly improved compared to a phase space grid approach which would require a complexity of (at least) O(N dx x N dv v ) for each time step where d x (resp. d v ) is the dimension in the spatial (resp. velocity) direction. Obviously, if an AP scheme is used, a similar time step would be considered between the grid and the particle approach. If we look at the comparison between PIC-BGK and MiMa, we observe in the following tabular that the cost of MiMa for one time step is about twice the cost of PIC-BGK and the complexity of both methods is O(N part ). Even if the complexity of the two methods is the same, PIC-BGK method is much more noisy in all the regimes. To reduce the noise, more particles are necessary compared to MiMa. In particular, when ε is small, PIC-BGK needs an important number of particles whereas MiMa gives very accurate results even if N part is very small (see Figure 16). For a given accuracy, MiMa is then faster than PIC-BGK. The CPU times are given in the following tabular to obtain results of Figure 18 for the electric energy. Note that the PIC-BGK simulation with ε = 10 −4 has been performed using an implicit Euler scheme for the equation on ω k , which enables to use a similar ∆t for MiMa and PIC-BGK. However, in this regime, MiMa needs very few particles so that the computational time is again in favor of MiMa. 6. Conclusion. In this work, an Asymptotic Preserving scheme using particles is proposed for the Vlasov-Poisson-BGK model in the fluid limit. This numerical scheme is based on a micro-macro decomposition; the AP property is ensured by using a suitable semi-implicit scheme to deal with the stiff source term. The main interest of the present approach compared to previous works on micro-macro decomposition is the fact that we use a particle method to discretize the micro part whereas a finite volume method is used to deal with the macro part. In this way, the numerical cost of simulations in the fluid regime is reduced. Moreover, in the non collisional limit, such an approach has the important property of reducing the numerical noise which is usually observed in the standard PIC approaches. This is a consequence of the micro-macro decomposition strategy and of the additional step which is performed to maintain this micro-macro structure along time.
Extensions of this approach to higher dimensions of the phase space are possible and will be the subject of future works. More complex collision operators (as the Landau operator of plasma physics) will also be considered combining this approach with relaxation techniques as in [22] for instance.
Appendix B. Computations relative to the second test case. We consider now the following test case: We have to compute an initial condition for the macro and for the micro equations.

ANAÏS CRESTETTO AND NICOLAS CROUSEILLES AND MOHAMMED LEMOU
Hence, the micro equation is initialized with: