Asymptotic preserving and time diminishing schemes for rarefied gas dynamic

. In this work, we introduce a new class of numerical schemes for rareﬁed gas dynamic problems described by collisional kinetic equations. The idea consists in reformulating the problem using a micro-macro decomposition and successively in solving the microscopic part by using asymptotic preserving Monte Carlo methods. We consider two types of decompositions, the ﬁrst leading to the Euler system of gas dynamics while the second to the Navier-Stokes equations for the macroscopic part. In addition, the particle method which solves the microscopic part is designed in such a way that the global scheme becomes computationally less expensive as the solution approaches the equilibrium state as opposite to standard methods for kinetic equations which computational cost increases with the number of interactions. At the same time, the statistical error due to the particle part of the solution decreases as the system approach the equilibrium state. This causes the method to degenerate to the sole solution of the macroscopic hydrodynamic equations (Euler or Navier-Stokes) in the limit of inﬁnite number of collisions. In a last part, we will show the behaviors of this new approach in comparisons to standard Monte Carlo techniques for solving the kinetic equation by testing it on diﬀerent problems which typically arise in rareﬁed gas dynamic simulations.

Abstract. In this work, we introduce a new class of numerical schemes for rarefied gas dynamic problems described by collisional kinetic equations. The idea consists in reformulating the problem using a micro-macro decomposition and successively in solving the microscopic part by using asymptotic preserving Monte Carlo methods. We consider two types of decompositions, the first leading to the Euler system of gas dynamics while the second to the Navier-Stokes equations for the macroscopic part. In addition, the particle method which solves the microscopic part is designed in such a way that the global scheme becomes computationally less expensive as the solution approaches the equilibrium state as opposite to standard methods for kinetic equations which computational cost increases with the number of interactions. At the same time, the statistical error due to the particle part of the solution decreases as the system approach the equilibrium state. This causes the method to degenerate to the sole solution of the macroscopic hydrodynamic equations (Euler or Navier-Stokes) in the limit of infinite number of collisions. In a last part, we will show the behaviors of this new approach in comparisons to standard Monte Carlo techniques for solving the kinetic equation by testing it on different problems which typically arise in rarefied gas dynamic simulations.
1. Introduction. The numerical simulation of kinetic equations involving many different applications ranging from rarefied gas dynamic and plasma physics to socio-economic models is a very active field of research. In this work, we focus our attention on the development of a new class of methods for rarefied gas dynamic problems. To this aim, concerning gas flow simulations, because of the intrinsic multiscale nature of many problems, it may happen that Euler or Navier-Stokes models are not sufficient in describing several different phenomena while kinetic models are very often the most adapted approaches. One of the key difference between fluid and kinetic models [10] is the high dimensionality of the mesoscopic approach which describes the state of the system by studying the time evolution of the probability of a particle to be in a given state in the six dimensional phase space at a given time [10,21]. On the other hand, the fluid models evolve macroscopic quantities such as density, temperature and mean velocity, which depend only on time and on the three dimensional physical space variables. It is a matter of fact that the kinetic description is much richer than the hydrodynamic one but the price to pay is very often too expensive numerical simulations which avoid the use of these approaches in practice.
For this reason, the numerical techniques for solving kinetic equations are often designed in such a way that the computational cost is as low as possible. However, it is undeniable that, even if many progresses have been done in the recent past [21], the goal of simulating realistic problems by means of kinetic models with deterministic techniques has not still be achieved. Thus, in practice, the most frequently used methods are based on probabilistic techniques, the most known method belonging to this category being the Direct Simulation Monte Carlo method (DSMC) [1,3,8,33,34]. However, this kind of approaches which work very well for stationary problems (thanks typically to time averaging techniques) are only poorly accurate if few particles are used or computationally too expensive if many particles are employed for unsteady problems. To overcome problems related to the low convergence rate and numerical noise of Monte Carlo schemes, several methods have been designed in the recent past. We quote in particular the review papers [8,35] for an overview on efficient and low variance Monte Carlo methods and some recent papers concerning the use of these techniques for rarefied gas dynamic problems [13,16,17,19,18,24,25,37,38] and plasma physics [11,6,7,9].
A second important drawback of kinetic approaches is that the collision term becomes stiff when interactions become important (i.e. when the system is close to the hydrodynamic limit). Thus, it turns out that standard explicit schemes lose their efficiency due to the necessity of using very small time steps to solve the collisional scale. In addition, in many situations, the fluid limit may occur in some regions of the domain and at some given times, while the kinetic regime is the most probable one in the rest of the domain. In these cases, one is normally obliged to resolve the micro scales in order to remain stable and consistent, but this requires very small time steps and phase space cells. On the other side, simulations have to be performed on macroscopic lengths, which makes the problem very challenging. Hence, domain decomposition approaches [5,14,15,31,27,28,22,36,23] can be adopted. However, the connection of the different models demands specific development as well as the interface identification is not always a simple task to solve. Thus, as an alternative, the problem has been recently addressed by designing the so-called asymptotic preserving schemes [4,12,20,26,29,30]. These methods are able to overcome the above numerical restrictions and automatically degenerate to consistent discretizations of the limiting models when the parameters which characterize the microscopic behaviors goes to zero. However, while the computational cost related to the solution of the microscale dynamics is in this way overcame, the dimensionality of the problem remains unchanged even in the limit case.
In this paper, we propose a new approach to solve kinetic collisional equations by using hybrid Monte Carlo techniques which addresses both the problem of the high AN APTD SCHEME FOR THE NAVIER-STOKES EQUATIONS 3 dimensionality and of large numerical noise. This work, inspired by [17,19,13], modifies a recent approach [11] in order to design a numerical scheme which is more efficient. In particular, as opposite to standard techniques, this new approach reduces the computational cost as the equilibrium state (which can be described by hydrodynamic like systems) is reached. In addition, the statistical error due to the particle part of the solution decreases as the number of interactions increases, realizing a variance reduction method which effectiveness depends on the regime studied.
The solution we propose is the following: we couple a Monte Carlo method for the solution of the kinetic equation with a finite volume method for the solution of the macroscopic equations in each point of the computational domain. Successively, we construct our algorithm in such a way that the computation of the microscopic solution can be avoided by using a simple asymptotic preserving (AP) method ; in such a way, the computational cost of the solution only depends at each time step on the quantity of solution which is out of the equilibrium state. In order to provide the correct solution for all the regimes, the macroscopic moment equations are coupled to the kinetic equation through a kinetic correction term, which takes into account departures from thermodynamical equilibrium. We both consider as a macroscopic equations the compressible Euler and the Navier-Stokes equations. As the equilibrium is approached, the number of particle in the Monte Carlo method is reduced which causes the computational cost to diminish and to be equivalent to the computational cost of a classic numerical method for the hydrodynamic equations in the limit. Thus, the method is automatically costly diminishing without imposing any artificial transition to pass from the microscopic to the macroscopic model at the contrary to domain decomposition techniques in which a transition region should be artificially imposed. In this sense, the schemes proposed here realize an automatic transition from kinetic to hydrodynamic which only depends on the real physics and not on numerical artifacts. Moreover, at the contrary to standard domain decomposition methods in which the cost is low only in regions of pure equilibrium, with these schemes the computational cost continuously diminishes passing from one regime to the other, while at the contrary to standard AP schemes the reduction of the complexity is not only due to the overcoming of the stiffness of the collisional scale, but also to the reduction of the dimensionality.
The remainder of the paper is organized as follows. In the next section we recall some basics about the collisional kinetic equations and their fluid limit. A reformulation of the kinetic equation which permits to design the numerical scheme is described in section 3, while the scheme itself is described in Section 4. Numerical results are presented in Section 5 which show the better performances of this approach with respect to classical Monte Carlo methods. Section 6 is used to draw some conclusions and suggest future developments. In Appendix are reported some detailed computations used through the paper.
2. The kinetic BGK equation, the compressible Euler and Navier-Stokes limits. We consider a simplified kinetic equation in which the interactions between particles are described by relaxation towards the local Maxwellian equilibrium [10]. This is the so-called BGK model and it is considered to be a valid alternative to the more complex Boltzmann operator in fluids which are not far from the thermodynamical equilibrium. It reads [10,21] where where ν = ν(ρ, T ) > 0 is a given relaxation frequency and measures the average time between two collisions, while ρ and T are the density and temperature of the gas defined below. We indicate the so-called local Maxwellian by M [f ] or M [U ] to stress the fact that it depends on the distribution function f but only through its moments U . It takes the form where u is the mean velocity. The vector U = (ρ, u, T ) is the vector of the macroscopic quantities obtained by integrating the distribution function over the velocity space multiplied by the so-called collision invariants where mf := R d f (v)m(v)dv and E = 1 2 ρu 2 + d 2 ρT is the total energy. Hence, M (U ) and f share the same first three moments in v. The kinetic equation is completed by boundary and initial conditions for f . When the gas is dense and temperature is sufficiently low, the Knudsen number is typically very small. In this case, the gas appears macroscopically in equilibrium. This is the fluid limit model [2,10] obtained taking the limit ε → 0 in (1). Let us briefly describe this limit. Multiplying (1) by m(v) and integrating with respect to v yields This is equivalent to the following system When ε goes to zero in (1), the distribution function tends to the local Maxwellian M (U ) given by (3). The previous system can then be closed, and the momentum (v ⊗ v)f = (v ⊗ v)M (U ) and energy fluxes |v| 2 vf = |v| 2 vM (U ) can be expressed as a function of U with p = ρT the pressure. Using the Chapman-Enskog expansion f = ∞ n=0 ε n f n , where the f n may depend on ε, higher order fluid models can be derived. We AN APTD SCHEME FOR THE NAVIER-STOKES EQUATIONS 5 specifically consider in this work the first order approximation to the distribution function f , which means we take f = f 0 + εf 1 . This choice gives after simple computations f 0 = M (U ) and f 1 = −(I − Π M )v ·∇ x M (U ), with Π M the orthogonal projection to be precised in the next Section. This leads to the usual compressible Navier-Stokes equations is the stress tensor and Q = k∇ x T the heat flux.
3. Derivation of micro-macro models. This section is devoted to the derivation of two different micro-macro models starting from equation (1) which are the basis of the numerical methods developed next. The first micro-macro model consider the compressible Euler equation as the macroscopic model while the second micro-macro decomposition employs the Navier-Stokes equations as a macroscopic model. Concerning the first case, the first step we write the distribution function f according to the following decomposition Now, since M = M (U ) and f shares the same first three moments, we have (recalling which in particular implies mg = 0. Let now T be the transport operator T f = v · ∇ x f , then the kinetic equation (1) writes Denoting now 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 (Q) the null space of the operator Q, the explicit expression of the projection operator can be analytically computed. This is given by (see [4]), for all ϕ ∈ L 2 (M −1 dv) the micro-macro model for the unknowns (g, U ), equivalent to the kinetic equation (1), can be written as follows Let us consider now the second decomposition. It reads Now, since M and f shares the same first three moments, we still have U (t, x) = mf = mM , which means mf 1 = mg = 0. Injecting the decomposition (10) into the kinetic equation (1) and applying the projection operator Π M gives which is equivalent to the following equation on the moments U where DU = ∇ x · vm(I − Π M )T M corresponds to the Navier-Stokes terms (right hand side of (5)). Now, injecting the decomposition (10) into (1) and applying (I − Π M ) leads to the microscopic part and using the definition of f 1 , we finally get the following model In the following section we introduce the numerical schemes based on the two cases: system (9) will be referred as Case 1 whereas system (11) will be referred as Case 2.

4.
A new Time Diminishing Asymptotic Preserving class of methods for kinetic equations. The goal is to design a class of schemes for solving kinetictype equations which avoids the resolution of the small scale dynamics induced by the particle interactions and which cost diminishes as the equilibrium state is approached. The class of scheme derived in this paper is based on a particle approach and on the following requests • It should be faster than standard deterministic methods (like finite volume or semi-Lagrangian) designed for solving kinetic equations. • The statistical error should be smaller than the one of standard Monte Carlo schemes. • The collisional scale should not impose time steps limitations.
• The computational cost should diminish as the equilibrium state is approached.
At the same time, the variance should diminish as the number of interactions increases. • The computational cost should be less than that of domain decomposition approaches. • The scheme should not need the introduction of artificial interfaces or criteria to pass from one regime to the other one.
To compute the solution of the kinetic equation for the perturbation g, a splitting procedure is performed (for both decomposition cases (9) or (11)) • Solve the transport part: ∂ t g + T g = 0.
• Solve the source part for the Case 1 (system (9)) or for Case 2 (system (11)), In both cases, during the second step of the above procedure particles are eliminated or respectively created depending on the shape of the solution at time n. If the amplitude of the perturbation g in one given spatial cell increases, the number of particles increases in order to better describe the departures from equilibrium. If, instead, the amplitude of the perturbation g in one given cell decreases, particles are eliminated. Thus the computational cost depends at each time step on the absolute value of the perturbation function g, no artificial conditions are imposed on the number of particles used, only the discretization of the source term determines the number of samples which are needed to describe the problem. In the following we details the algorithm described for both cases.

Solution of the kinetic equation for the perturbation function.
We consider the solution of the equations (9a) and (11a) in a time interval [0, ∆t] by using Monte Carlo methods. This approximated solution is given by an operator splitting between transport ∂ t g + T g = 0, (12) and the source terms. They read respectively for the first and the second cases In general, high order time splitting schemes can also be used for Monte Carlo approaches. However, in this work we did not consider this possibility. In the sequel, we restrict to the one-dimensional case in x and v. Let us now introduce a space discretization of the interval [x max , x min ]: We also consider a time discretization of time step ∆t: t n = n∆t, n ≥ 0. The discretization of the space domain is not needed for the transport step which is solved exactly by simply pushing forward in time the particles. However, it is necessary to solve the source part which acts locally in space by reconstruction of the moments of the perturbation g in the space cells. In a Monte Carlo method, the distribution function g is approximated by a finite set of N particles where x k (t) represents the position, v k (t) the velocity and ω k the weight of each particle. Note that for the particles positions and velocities, we use the notation x k (t) and v k (t) whereas x j and v j are used for the fixed phase space mesh. The initial position, the initial velocity and the weight are defined as following. Let N p be the total number of particles used to discretize the original problem (1), where the unknown is the distribution function f . Then the mass m p of each particle is simply defined as In the same way, we define N 0 j the number of initial particles in each cell I j (centered around x j ) which is given by in particular we have Nx j=1 N 0 j = N. Then, one computes the number of particles with positive weight ω k = m p and negative weight ω k = −m p by computing the integer numbers (17) and (18) should be interpreted through a stochastic rounding procedure defined as with ⌊x⌋ the integer part of x. Once the initial number of samples N 0 j is fixed, positions (x k ) are uniformly assigned over the cells centered in x j . Successively, fixing a grid in velocity space of size ∆v and the boundaries of the velocity space v min and v max , we approximate g(t = 0, x j , ·) by a piecewise constant function on each velocity cell of size ∆v. Then, in each cell x j , we have to sample N 0 j particles in the velocity direction. To do so, the velocities v k (0) are assigned by inverting the following equation The particle velocities v k have a positive weight ω k = m p when they are extracted from G + , and a negative weight ω k = −m p when they are extracted from G − . The reals ξ k are random numbers in the interval We discuss now the two steps of the splitting. Suppose that g n , M n and f n 1 are known as well as x n k and v n k for k = 1, . . . , N n g , where N n g = Nx j=0 N n j and During the transport step, we solve the characteristic equations between t n and t n+1 which means that particles are pushed forward in time. Any type of time integrator can be employed to solve this step, we consider the standard first order Euler solver which, knowing x n k , v n k reads x * k = x n k + ∆t v n k . Obviously, more accurate solvers can be used, we refer to [21] for an overview of time integration techniques. Before going to the relaxation part, we perform an intermediate step in order to reduce the number of particles in each cell as much as possible. This intermediate step is as follows. Let N * j,v k the number of particles in the cell j having velocity v k after the transport part, and let N * ,± j,v k the corresponding number of particles having velocity v k with mass ±m p . Among this set of N * j,v k particles, we keep only since particles with opposite masses and same velocity can be eliminated. This is an important step since otherwise the number of particles may grow indefinitely in time especially for large values of ε, without providing any additional information on the distribution function. Now, we deal with the second part of the splitting. This is done by using an implicit-explicit scheme which permits to deal with the stiff terms. In the first case (9), it reads which gives In the second case (11), it reads which gives In the above equations, g * represents the solution after the transport step and the precise expression for are detailed in Appendix A. Equations (21)- (22) are convex combinations of two terms for each choice of the time step ∆t. The first term g * is a sink, the second one P or P 1 is a source. From a Monte Carlo point of view, these equations state that with probability (1 − q) = ∆t ε/ν+∆t , particles are discarded and with a probability q = ε/ν ε/ν+∆t , particles are created from a sampling of the sources P(x, v) or P 1 (x, v). More precisely, after the transport step, the functions P(x j , v) and P 1 (x j , v) are computed in each cell I j . After the pushing step (20), we first need to evaluate the cell to which each particle belongs to. Then, the algorithm to approximate (21) or (22) is the following by sampling P + (x j , v) for positive weights particles and P − (x j , v) for negative weight particles. This gives the total new number of particles in the cell The second case is considered by simply changing P by P 1 given by (24). As for the function g, we defined P ± (x j , v) by P ± (x j , v) = 1 2 [P ± (x j , v) ± P ± (x j , v)] (and similarly for P ± 1 (x j , v)). The positions for particles which do not disappear remain unchanged x n+1 k = x * k while the new positions x n+1 k for the new samples are computed uniformly over the space cell. Concerning the velocities, they remain the same if the particle is kept v n+1 k = v * k while their values depend on P(x, v) and P 1 (x, v) in the other case.
In the asymptotic regime ε → 0, we have q = O(ε) and M ± j = O(ε) since P(x, v) = O(ε) and P 1 (x, v) = O(ε) thanks to (23) and (24). As a consequence, the number of particles in the cell I j used to sample g satisfies which means that the number of particles diminishes automatically as ε → 0.
Note that a similar approach has been developed in [37] but only in the framework of Case 1 where the micro part of the micro-macro decomposition is approximated by an hybrid Monte Carlo type method with negative particles (see [38]).

4.2.
Numerical discretization of the macroscopic equations. In this section we discuss the discretization of the moment equations (9b) The first system is composed by an equilibrium part which corresponds to the classical compressible Euler equations and by a non equilibrium part which takes into account departures from equilibrium. The second system is composed by the same equilibrium part which gives the compressible Euler fluxes and by a diffusion part which gives the Navier-Stokes fluxes. The time discretizations are the following and For the space discretization of the Euler fluxes F (U ) in both systems, we consider second order MUSCL central schemes. They read with ϕ a modified slope limiter, α equal to the larger eigenvalue of the Euler system and χ n,± j the ratio of consecutive discrete gradients of the fluxes [32]. The slope limiters proposed here intend to avoid the increase of the statistical fluctuations which may occur by using high order reconstructions, in particular when the perturbation g is large. In such a case, it is better to smooth out the fluctuations by passing from a second order scheme to a first order one. Thus, the limiters are where ϕ c (χ n,± j ) is the usual slope limiter (we used the Van Leer one), and β(χ n,± j ) is the unit step function which depends onχ n,± j , the ratio of the discrete gradients of the equilibrium and non equilibrium fluxes. In details, when the non equilibrium fluxes exceed the compressible Euler fluxes of a certain size, then, automatically, the discrete second order fluxes pass to discrete first order ones. In the numerical test section, we fixed this threshold to 0.25.
We now discuss how to discretize the non equilibrium terms ∇ x · vmg n and the Navier-Stokes term εDU n . The Navier-Stokes term is discretized by employing the same second order MUSCL scheme used for the Euler fluxes where the numerical diffusion α is taken equal to zero. Concerning the perturbation terms, we introduce a filter to eliminate eventual fluctuations. These fluctuations are due to the moments a j = vg j and b j = v 2 g j which are not exactly zero on each cell j. To make the values of these moments as small as possible, we choose the weighted moving average method, which consists of a convolution of the pointwise values of these quantities a j and b j with a fixed weighting function, that is we replace a j and b j bȳ Once the filter imposed, a first order discretization is employed for the non equilibrium term without numerical diffusion, i.e. α = 0. We then define the numerical fluxes as follows 4.3. Projection step. Let us first notice that the exact preservation of the moments mg is not guaranteed by the scheme described above. Indeed, due to the use of a Monte Carlo method, even if the moments of g are equal to 0 at the beginning of the computation, i.e. for t = 0, we can not ensure that the property mg = 0 propagates in time. Let observe that this is due to the stochastic component of the solution. In fact, in the case in which a finite volume method is used for discretizing the kinetic equations (9) or (11) and employing the same time discretization used in the previous paragraph, the condition mg = 0 is fulfilled for all times if it is satisfied at t = 0. For this reason, we discuss in this section a way to project the function g in order to ensure that the condition mg = 0 is exactly fulfilled for all time in the case in which particles are used to approach the distribution g. It has been observed in [11] that this projection step allows a significant reduction of the statistical noise in the context of plasmas. However, in the present context, we noticed from numerical experiments that the enforcement of this property at a discrete level does not improve the numerical results while it may induce an additional computational cost.
Let us describe the projection step in our context. This is done locally on the cells and we therefore only consider one fixed spatial cell, the same procedure being repeated for all cells. In order to proceed, we suppose that the solution g n is known everywhere thanks to the knowledge of the distribution of particles with weight ω k , velocity v k and position x k , for k = 1, . . . , N g (with N g = Nx j=0 N j and N j = N + j + N − j , forgetting the time dependency), after using for one time step the scheme presented in section 4.1 from t n−1 to t n . In this setting, the property that we want ensure is Nj k=1 m(v k )ω k = 0 (31) in each cell centered around x j . In the sequel, we discuss how to restore the three conservations.

4.3.1.
Conservation of g . The loss of conservation is due to the first part of the splitting scheme for the kinetic equation, namely (12). In the sequel, we detail a way to ensure that the zero order moment of g is ensured: g = 0 in each spatial cell for all times.
If the density of g particles is not zero, this means that Nj k=1 ω k = 0 in one given cell centered around x j . This can be also interpreted by saying that N + j = N − j : the number of samples with positive weight is different from N − j , the number of samples with negative weight. Thus, in order to restore equilibrium between positive and negative weight samples, we can choose to eliminate some samples or to create new ones. The best choice is the following: if N + j > N − j we eliminate (N + j − N − j )/2 samples with positive weight and we create (N + j − N − j )/2 samples with negative weight. If conversely N − j > N + j we eliminate (N − j − N + j )/2 samples with negative weight and we create (N − j −N + j )/2 samples with positive weight. The elimination is done uniformly while for a creation of a new particle we proceed as following. First, we search from the discarded particles obtained from the solution of equation (21) if we have enough particles to restore the equilibrium between positive and negative weights. In this case, samples are taken with uniform probability from the positive or negative discarded samples. If conversely there are not enough samples from the discarded set computed in (21) or (22), new samples are created by replicating existing ones. When a new particle is created, its position is chosen uniformly in the cell grid. Finally, let observe that |N − j + N + j | can be an odd number, in this case the last particle we need to restore equilibrium is taken with the same probability from the positive or negative weight sets.

4.3.2.
Conservation of the first and second order moment of g. As for the conservation of the difference between positive and negative particle in one cell for all times, the conservation of the momentum and energy of the distribution g is due to the use of a particle approach and to the transport step (12) which modifies the position of the particles in the domain. In order to ensure conservation of momentum and energy, i.e. vg = 0, |v| 2 g = 0 we propose a modification of the renormalization used in [14]. The strategy is composed of two steps: in the first one we adjust the average momentum by simply translating the samples in the velocity space for each direction. In the second part, we adjust the energy by using two different scalings for the positive and negative samples. In details, the transformation is the following. Suppose to be known v n+1 k , which is recalled in this paragraphṽ n+1 k , the velocities after the transport and the relaxation steps described in (21). Suppose also that the mass conservation process described in the previous paragraph in each spatial cell has been done. Then the momentum conservation is ensured thanks to the following transformation v n+1 which is done for all particles which belongs to the cell centered around x j . Then, after computing the energy associated with the positive and the negative samples as the following transformation is imposed where c + and c − are defined by This is enough to ensure that the energy associated to the perturbation g is zero in each cell at each instant of time. We conclude this section by briefly describing the complete algorithm.

4.4.
Algorithm. We first give the algorithm in the case of decomposition (9). Suppose that two grids of width ∆x and a time step ∆t have been fixed. Suppose also that the moments U n and the functions g n are known in each point of the spatial grid at time n as well as particle number N n g , positions x k , velocities v k and weights ω k . Then the algorithm for passing from time t n to time t n+1 is the following: 1. Push forward the particles by solving (20). Compute ∇ x · vmg using (29) and (30). 8. Advance the macroscopic moments equation using (25). Let us now consider the case of decomposition (11). As before, let us suppose two grids of width ∆x and a time step ∆t have been fixed. Suppose also that the moments U n and the functions f n 1 are known in each cell of the spatial grid at time n as well as particle number N n g , positions x k , velocities v k and weights ω k . Then the algorithm for the decomposition (11) for passing from time t n to time t n+1 is the following 1. Push forward the particles by solving (20). 5. Numerical results. In this section we discuss several numerical results obtained with the schemes presented. From now on, the Asymptotic Preserving and Time Diminishing methods described in the Algorithm 4.4 will be called respectively APTD for the decomposition (9) and APTDNS (Navier-Stokes) for decomposition (11). The two approaches will be compared with a classical Monte Carlo method for the kinetic BGK equation (called MC). In order to have a reference solution we also use a deterministic second order in space scheme based on a Discrete Velocity Model formulation (DVM scheme) of equation (1) and a finite volume method for the hydrodynamic model (Euler scheme). For all tests the number of cells in the physical space will be the same for the five cited methods, while for the APTD, the APTDNS and the MC methods the mass of a test particle will be the same. As clear from the description of the APTD and APTDNS methods in the previous Sections, given the mass of a particle for the MC method m p , the mass of a particle for the APTD or the APTDNS methods will be either m p or −m p depending on the sign of g. The Section is divided into three parts dedicated to test our method against a smooth solution and against some classical problems arising in gas dynamics. 5.1. Test 1: Unsteady Shock. In this test, a gas with constant density (ρ = 1) and temperature T = 5 is pushed against a wall at a speed u = −1 with perfectly reflecting boundary conditions. In the fluid limit such an initialization produces a shock wave which moves from the left to the right of the domain, after the shock, the speed of the flow being equal to zero. When collisions are not sufficient to reach the equilibrium state at each instant of time, the shock wave is smoothed out and a continuous transition from the right to the left state of the macroscopic quantities is realized which is faster as the fluid limit is approached. The number of cells in physical space is equal to N x = 300, the final time is set equal to t f = 0.067 while the test has been run for ε = 10 −2 , 10 −3 , 10 −4 . The number of particles used in the Monte Carlo method is at the beginning 10 6 and it grows with time since the density is physically non conserved in the domain for this test case, instead it grows due to the formation of the shock wave. For the APTD and APTDNS methods, the mass of each particle is the same of the mass of the particles in the original MC method (m p ≃ 10 −7 ). Concerning the DVM and the Euler schemes, a second order MUSCL scheme has been used in both cases with a number of points in velocity space for the DVM equal to 120.
In Figure 1, the density, the mean velocity and the temperature for respectively the APTDNS, on the left and the MC scheme, on the right, are reported for a rarefied regime which corresponds to ε = 10 −2 . In Figure 2, the same macroscopic quantities are reported for ε = 10 −3 which corresponds to the intermediate regime. Finally in Figure 3, these quantites are reported for a regime close to the fluid limit, ε = 10 −4 . The results show as expected that the reduction of the variance is more important when the interactions between the molecules become more important. For the three different interaction scales chosen, an important reduction of the statistical noise with respect a classical Monte Carlo scheme is observed. In Figure  4 (left), we can observe that the number of particles effectively employed in the APTD method in comparison with the number of particles used by Monte Carlo. On the middle part of Figure 4, we report the ratio between these two numbers where we clearly see that as expected this number decreases as the equilibrium is approached. In particular we observe that, in the worst studied scenario (ε = 10 −2 ), the number of particles effectively used is less than 1.5% of the number of samples of the corresponding MC method, while this ratio goes to 0.01% for ε = 10 −4 . This gives precise indications concerning the computational cost of the method, since most of the computations are due to the kinetic part of the model, being the computations due to the macroscopic part of the solution negligible. In Figure 4 (right), we finally report the ratio between the number of particles used by the  APTD and the APTDNS schemes as a function of time for the different values of the Knudsen number considered. We see that this ratio as expected is always less than 1 and in particular we see that in the case of small Knudsen number, the number of particles used in the APTDNS is around 25% less than the number of particles of the APTD scheme.    For the APTD and APTDNS methods the mass of each particle is the same of the mass of the particles in the original MC method, which makes m p ≃ 10 −7 . Concerning the DVM and the Euler schemes a second order MUSCL scheme has been used in both cases with a number of points in velocity space for the DVM equal to 120.
In Figures 5, 6 and 7, the density, the mean velocity and the temperature for respectively the APTDNS, on the left and the MC scheme, on the right, are reported for the same regimes considered in the first test: ε = 10 −2 , 10 −3 , ε = 10 −4 . The results show as before that the reduction of the statistical error is more important when the interactions grow. For all regimes considered, our schemes give less noisy results than Monte Carlo. Finally, Figure 8 (left part) shows the number of particles used by the APTD method in comparison with the number of particles used by the Monte Carlo one. On the middle of Figure 8, the ratio between these two numbers is reported while on the right part, the ratio between the number of particles used by the APTD and the APTDNS schemes is shown. The results can be resumed by saying that the APTD methods produce less noisy solutions with less particles and this number decreases when the solution becomes closer to the fluid regime. In particular, the APTD method which employs the second decomposition is able to produce solutions with less particles than the APTD method since the Navier-Stokes terms are computed in a deterministic way.  this reason, we do not report the results for this test. Instead, in Figure 9, we report the L 1 norm of the error for APTD (left column) and MC (right column) for the different Knudsen numbers. These have been obtained by measuring the differences between the macroscopic quantities computed by APTDNS and MC schemes with respect to a reference solution. The reference solution has been obtained by a second order in space deterministic scheme based on discrete velocity model version of the original kinetic equations. The results show that we have a reduction of the error compared to classical MC schemes for all the situations considered. In particular, the reduction of the error becomes larger as the fluid limit is approached. 6. Conclusion. In this work, we have presented a new numerical method for the BGK equation which enjoys the following properties: (i) its statistical noise is smaller than the one of standard Monte Carlo methods; (ii) it is asymptotically stable with respect to the Knudsen number; (iii) its computational cost as well as its variance diminish as the equilibrium is approached; (iv) no artificial criteria is required. The method is based on a micro-macro decomposition (Euler-kinetic or Navier-Stokes-kinetic) for which the macro part is solved using a finite volume method whereas the micro part uses a Monte-Carlo method. This enables to derive a low variance Monte Carlo based method for which additionally, the number of particles used to sample the micro unknown diminishes automatically as the fluid regime is approached. The numerical results illustrate the efficiency of the proposed method compared to the standard Monte Carlo approach.
In the future, we would like to extend the present approach to the case of the Boltzmann operator, to the solution of the Vlasov equation and to the challenging case of diffusive limits.