The Moment Guided Monte Carlo method for the Boltzmann equation

In this work we propose a generalization of the Moment Guided Monte Carlo method developed in [11]. This approach permits to reduce the variance of the particle methods through a matching with a set of suitable macroscopic moment equations. In order to guarantee that the moment equations provide the correct solutions, they are coupled to the kinetic equation through a non equilibrium term. Here, at the contrary to the previous work in which we considered the simplified BGK operator, we deal with the full Boltzmann operator. Moreover, we introduce an hybrid setting which permits to entirely remove the resolution of the kinetic equation in the limit of infinite number of collisions and to consider only the solution of the compressible Euler equation. This modification additionally reduce the statistical error with respect to our previous work and permits to perform simulations of non equilibrium gases using only a few number of particles. We show at the end of the paper several numerical tests which prove the efficiency and the low level of numerical noise of the method.


Introduction
The kinetic description of a gas is based on the integro-differential Boltzmann equation [2,8]. This equation is satisfied by the probability distribution function of the particles depending on the time t ≥ 0, the space variable x ∈ R d and the velocity v ∈ R d . In addition to the high dimensionality of the problem, the Boltzmann collision term that characterizes the kinetic equation is a five fold integral very hard to treat due to its nonlinear nature and to the physical properties which need to be preserved during its resolution.
For these reasons, the numerical simulation of the Boltzmann equation with deterministic techniques appears to be very difficult and thus, in practice, probabilistic techniques such as Direct Simulation Monte Carlo (DSMC) methods are largely used in real situations due to their flexibility and low computational cost compared to finite volume or spectral methods for kinetic equations [1,2,6,7,28,29]. However, DSMC solutions are affected by large fluctuations. Moreover, in non stationary situations the impossibility to use time averages to reduce fluctuations leads to poorly accurate solutions or computationally very expensive simulations. To overcome this problem, several methods have been developed. We quote [6] for an overview on both efficient and low variance Monte Carlo methods. For applications of variance reduction techniques to kinetic equation we mention the works of Homolle and Hadjiconstantinou [22] and [23]. We mention also the research papers by Pullin [30] which developed a low diffusion particle method for simulating compressible inviscid flows and the work of Weinan E. and co-authors [20] on the development of multiscale numerical methods. We quote also the work of Lemou and co-authors who developed a low variance method for the Vlasov equation close to the fluid limit [9]. We finally recall the results obtained by the author of the present paper [15,16] on the construction of efficient and low variance methods for kinetic equations in transitional regimes.
A second important drawback of kinetic approaches is that the collision term becomes stiff close to the fluid regimes. A non dimensional measure of the importance of collision is given by the Knudsen number which is large in the rarefied regions and small in the fluid ones. Thus, standard computational approaches lose their efficiency due to the necessity of using very small time steps in deterministic schemes or, equivalently, a large number of collisions in probabilistic approaches. One possibility which permits to avoid the severe time step restrictions caused by the small free path limit is to construct the so-called Asymptotic Preserving schemes [9,18,24] which permits to avoid time steps dependence of the mean free path. Another possibility to avoid the excessive computational cost and the stiff regimes is to construct numerical schemes which combine continuum models with microscopic kinetic models [4,5,10,12,13,14,21,26,31,32]. However, this approach has several drawbacks, first the identification of the different regions is not trivial and second match the two models at the interfaces is not simple. For these reasons alternatively approaches are often preferred.
The two main goals of this work are first, design a particle closure of the macroscopic fluid equations which avoid the severe time step restrictions caused by the collisional scale and second to reduce the statistical error caused by the use of DSMC methods. This approach has explained in detail in [13] can also be the basis of very efficient domain decomposition strategies. In our first paper [11], we developed a method for the simple case of the BGK equation which consists in forcing particles to match prescribed sets of moments given by the solution of deterministic equations. Now, we generalize the approach to the case of the Boltzmann equation. Moreover, we construct the scheme in such a way that the time step restriction due to the small Knudsen number are avoided, thanks to the introduction of exponential schemes for the resolution of the Boltzmann collision term [17]. In particular, in the limit ε → 0 the scheme becomes automatically an high order method for the compressible Euler equation in which the numerical noise completely disappears. In this sense, the method belongs to the class of the so-called asymptotic preserving methods [17,18,19,24].
The method is based on the following idea. The resolution of the problem is done through the kinetic equation which is solved by DSMC method and by a set of closed moment equations which are solved by means of finite volume techniques. In order to provide the correct solution for all the regimes of the Knudsen number, the moment equations are coupled to the Boltzmann equation through a kinetic correction term, which takes into account departures from thermodynamical equilibrium. In this sense, we perform a closure to the infinite set of moment equations, using the solution of the kinetic equation. Observe that the kinetic model will be sufficient to give the correct solution to the problem. However, as already explained, the Boltzmann equation should be resolved by very fast methods because of its too high computational cost. Typically in real situations Monte Carlo is the only possible choice. The drawback is that the solution is only poorly accurate. On the other hand, the moment equations can be resolved by using high order finite volume techniques. The result is a method which has a computational cost comparable to the cost of DSMC method but which a smaller level of numerical noise. Moreover, the stiffness of the problem is removed thanks to the exponential resolution of the collision term. This gives a method which is faster than classical DSMC methods [1,2] in the limit of the Knudsen number which goes to zero i.e. in the fluid limit. Finally, in order to recover at each time step the same moments from the solution of the kinetic equation and from the solution of the moments equations and then advance in time, we constrain the DSMC method to match the moments obtained through the deterministic resolution of the macroscopic equations in such a way that the higher accuracy of the moments resolution improves the accuracy of the DSMC method. We experimentally show that this is indeed the case.
The remainder of the paper is organized as follows. In the next section we recall some basic notions on the Boltzmann equations and its fluid limit. The scheme is described in section 3. In particular, the numerical method used for the kinetic equation is described in section 3.1, the numerical method used for the moment equation is described in section 3.2, while the matching procedure is described in section 3.3. In section 4 numerical examples which demonstrate the capability of the method are presented. Finally some conclusions and remarks are drawn in the last section.

The Boltzmann equation and its fluid limit
We consider the Boltzmann equation of rarefied gas dynamics [8] Here f (x, v, t) is a non negative function describing the time evolution of the distribution of particles with velocity v ∈ R 3 and position x ∈ Ω ⊂ R dx at time t > 0. In the sequel for notation simplicity we will omit the dependence of f from the independent variables x, v, t unless strictly necessary.
The operator Q(f, f ) which describes particle interactions has the form where and B(|v − v * |, n) is a nonnegative collision kernel characterizing the microscopic details of the collision given by with γ ∈ [0, 3) and σ the collision cross section which in general depends on the collision angle. The case γ = 1 is referred to as hard spheres case, whereas the simplified situation γ = 0, is referred to as Maxwell case. Note that in most applications the angle dependence is ignored and σ is assumed constant [2].
The operator Q(f, f ) is such that the local conservation properties are satisfied where m(v) = (1, v, |v| 2 2 ) are the collision invariants. In addition it satisfies the entropy inequality Integrating (1) against its invariants in the velocity space leads to the following set of non closed conservations laws Equilibrium functions for the operator Q(f, f ) (i.e. solutions of Q(f, f ) = 0) are local Maxwellian of the form where ρ, u, T are the density, mean velocity and temperature of the gas in the x-position and at time t defined as We will denote by U = (ρ, ρu, E), When the number of collisions is large which means that the mean free path between particles is very small, it is convenient to rescale the space and time variables in (1) as which leads to where ε is the so-called Knudsen number proportional to the mean free path while the primes have been omitted for simplicity. If we pass formally to the limit for ε → 0 we get f → M [U ]. This leads to the well-known hyperbolic system of compressible Euler equations for the macroscopic variables U where I is the identity matrix.

The Monte Carlo method for the Boltzmann equation
The method is based on the decomposition In the above formula the function g represents the non-equilibrium part of the distribution function.
Since f and M [U ] share the same first three moments we get which means that the non-equilibrium part g has always non positive parts. With the above definitions U , f and g satisfy the following system of equations The proof of the above statement can be found in [12], where we also remind for details. As explained in the introduction, the goal of the method is to solve equation (18), which means to close the system of moments equations. In order to accomplish that task, we need to know the time evolution of g, which is given by the solution of equation (19). Of course solving equation (19) will be sufficient to compute the solution also of (18). However, the numerical solution of equation (19) is very expensive and thus we seek for very fast solvers, as for instance the Monte Carlo method. The counterpart of fast solvers is that the solution of (19) will be normally known with very low accuracy. Thus, our statement is the following: if the departure from equilibrium g is small, the low accuracy and the high level of numerical noise with which g is known from (19) will still permit to compute the moments in (18) with an error which will be smaller than the error with which equation (19) is solved. We will numerically show that this is indeed the case. In particular, thanks to the way in which the collisional operator is solved, the coupling method will permit to make disappear g from the computation of the moments in the limit ε → 0. The precise time marching procedure will be detailed next, the method can be summarized as following: at each time step t n 1. Solve the kinetic equation (19) with a Monte Carlo scheme and obtain the set of moments U * = mf * .
2. Solve the fluid equation (18) with a finite volume/difference scheme using the particles solution to close the system and obtain a second set of moments U n+1 . This means evaluate the term ∇ x · vmg with particles.
3. Match the moments of the kinetic solution with those of the fluid solution through a linear transformation of the samples velocity f n+1 = T (f * ) so that mf n+1 = U n+1 and continue.
In the first step, the method can be generalized by substituting the Monte Carlo method with any low accurate but very fast solver. The other two steps are the key points of the method they involve the evaluation of ∇ x · vmg and the matching procedure detailed next. We observe that it is possible to improve the method, adding to system (18) additional equations for the time evolution of higher order moments and get with m n = v n and n ≥ 3. However, this will lead to more complex equations which discretiazion will not be easy. We will not discuss this possibility here.

Solution of Boltzmann equation by exponential Runge-Kutta Monte Carlo methods
The starting point in the solution of kinetic equations by Monte Carlo methods is given by an operator splitting of (1) or equivalently (19) in a time interval [0, ∆t] between relaxation and free transport Even if this splitting is limited to first order it permits to treat separately the hyperbolic free transport from the stiff relaxation step. Observe that high order splitting schemes are possible with Monte Carlo techniques. However, we did not consider this possibility in this work, we remind to the future for the exploration of high order time discretizations.
We introduce now a space discretization of mesh size ∆x and a time discretization of time step ∆t. The discretization of the domain is not needed for the transport step which is solved exactly by pushing the particles. However, it is necessary to solve the collision part of the problem which acts locally in space and thus it turns to be necessary to solve the full problem. In Monte Carlo simulations the distribution function f is discretized by a finite set of N particles where X i (t) represents the particle position in the three spatial directions, V i (t) the particle velocities in the velocity space, and α i (t) the weight to associate to each particle. The constant m p is defined at the beginning of the computation in the following way During the transport step (23), the particles move to their next positions according to where (24) with (26) and a constant α i (t) represents an exact solution of equation (23). The role of the function α will be clarified in the moment matching procedure.
We consider now the solution of the collision step (22). The effect of this step is to change the shape of the velocity distribution and to project it towards the equilibrium distribution leaving unchanged mass, momentum and energy. The collision operator acts locally in space which means that we solve it independently in each spatial cell. The particles are assumed to have all the same weight in one cell. The initial data of this step is given by the solution of the transport step because of the splitting.
There exists many different ways to solve the collisions, see for instance [1,2,8], however, we seek for a method which will be stable for all choices of ∆t independently of ε. This will permit to avoid the stiffness of the equation and at the same time to reduce the statistical fluctuations. Since we aim at developing unconditionally stable schemes, the most natural choice would be to use implicit solvers applied to (22). Unfortunately the use of fully implicit schemes for (22) is unpracticable in the case of DSMC methods. To overcome these difficulties, we introduce a method based on exponential integrators. First we rewrite the homogeneous equation (22) in the form where P (f, f ) = Q(f, f ) + µf and µ > 0 is a constant such that P (f, f ) ≥ 0. Typically µ is an estimate of the largest rate of the negative term in the Boltzmann operator. For example for Maxwellian molecules we have and µ = ̺. By construction the following property holds this means that P (f, f )/µ is a density function. The homogeneous equation can be written now, adding and subtracting M [U ], in the form Note that even if M [U ] is nonlinear in f thanks to the conservation properties (6), it remains constant during the relaxation process.
We apply now to the reformulated equation (30) an approach based on the so-called exponential integrators where the exact solution of the linear part is used for the construction of the numerical method [17]. In order to derive the methods it is useful to rewrite (30) as The above form is readily obtained if one multiplies (30) by the integrating factor exp(µt/ε) and takes into account the fact that M [U ] does not depend of time during the collisional process. A class of explicit exponential Runge-Kutta schemes is then obtained by direct application of an explicit Runge-Kutta method to (31). More generally we can consider the family of methods characterized by where f * is the distribution function value after the transport step, F (i) are called stages, c i ≥ 0, while the coefficients A ij and the weights W i are such that with coefficients a ij , c i and weights w i given by a standard explicit Runge-Kutta method called the underlying method. Various schemes come from the different choices of the underlying method. The most popular approach is the integrating factor (IF) method. For the so-called Integrating Factor methods, which correspond to a direct application of the underlying method to (31), we have with λ = µ∆t/ε. In the sequel, among all possible choices for discretizing (27), we will use the first order in time integrating factor scheme which reads which is based on the simple firs order in time explicit Euler Runge-Kutta method. Observe that this method permits to avoid the stiffness of the collisional operator. This is, in fact, unconditionally stable for all choices of time step ∆t. At the same time, the method proposed is, in the limit ε → 0, nothing else than a kinetic scheme for the numerical solution of the compressible Euler equation.
Observe, in fact, that when ε → 0 at each time step the distribution function is projected on the The Monte Carlo interpretation of equation (35) is the following. If we define with probability A the velocity of the particle does not change during the collision step. With probability B the particle occurs in one collision, the details of the collision are described by the operator P , this is, for instance, the gain part of the collisional process, described by Bird [2]. With probability C the velocity of the particle is replaced by a particle sampled from the Maxwellian distribution M [U ]. We will discuss in the matching section the way in which we link the solution of the kinetic equation to the solution of the moment equations.

Solution of the moments equations
In this section we detail the way in which the moments equations have been discretized. We suppose for the construction of the scheme that the function g n is known at time n. The numerical scheme proposed will take advantage from the fact that part of these moment equations are in fact nothing else that the compressible Euler equations We first solve the set of compressible Euler equations and then we consider the discretization of the kinetic flux ∇ x · vmg . For the space discretization of the equilibrium fluxes we use a second order MUSCL central scheme. For simplicity, we indicate in the same way the numerical flux in one or in more spatial dimensions. Thus, we have The discrete flux reads [25,27] with ϕ the Van Leer slope limiter finally, the variable χ ± is defined as following where the above ratio of vectors is defined componentwise and α is equal to the largest eigenvalue of the Euler system. We now discuss how to discretize the non equilibrium term ∇ x · < vmg >. To this aim, the first step is to introduce a filter to eliminate some fluctuations. We choose the weighted moving average method, which is a convolution of the pointwise value of g with a fixed weighting function. Thus, given the function g n , we define The smoothing permits to reduce the fluctuations but on the other hand it causes a degradation of the accuracy in the solution. Thus, in practice, we use a weak filter in order to keep the solution as precise as possible. In the numerical simulations, we used K = 1, ω −1 = ω 1 = 1/6 and ω 0 = 2/3. Observe that a different smoothing can be used at the particle level instead that at the moments level. This different smoothing has to be applied during the reconstruction of the moments of g from the particles. We discuss this possibility in the next section. Once that the filer is applied, the non equilibrium term is discretized with the same second order MUSCL scheme used for computing the flux of the Euler equations where, however, the numerical diffusion is taken equal to zero. This is because this term is at the leading order a diffusion term, and thus it does not need numerical diffusion for stability to be assured. The complete scheme for the macroscopic equations reads where Ψ j+1/2 (< vm g n >) is defined as Ψ j+1/2 (< vm g n >) = 1 2 (< vm g n j > + < vmg n j+1 >) + 1 4 (σ n g,j − σ n g,j+1 ) with σ n g,j = < vm g n j+1 > − < vm g n j > ϕ(χ n g,j ) and where again the above ratio of vectors is defined componentwise.

The Moment Matching
In this section we discuss the details of the moment matching and the coupling between the kinetic and the macroscopic model. Suppose, known f n−1 , g n−1 and U n−1 . The time marching procedure is the following 1. Solve the moments equations from time n − 1 to time n using equations (44-45) this gives U n .
2. Solve the transport part of the kinetic equation from time n − 1 to the intermediate stage * pushing the particles as in equation (26).
3. Compute the moments of the Boltzmann equation after the transport of the particles. This gives U n , the collision part, being conservative, does not alter the moments. The reconstruction of the moments from the particles will be discussed next. 8. Plug the value of g n in the moments equations and compute U n+1 using (44-45) and continue.
We discuss now, how to reconstruct the moments from the particles (point 3 of the above procedure), the matching of the different moments (point 4) and the computations of g n from f n (point 6).
We first consider the matching of the mass. To this aim, let consider the set of particles X 1 , . . . , X NI j inside the cell I j after the transport step. The corresponding mass is computed as others possible reconstructions will be discussed next. In our previous paper, [11], among the possible techniques that can be used to restore a prescribed density we choose to replicate or discard particles inside the cells. Here, in order to restore the mass we assign to each particle inside the cell j the new value where ̺ n j is density computed from the solution of the moments equations. This rescaling of the mass introduce an error in the evaluation of the distribution function f . However, we recall that the two models, the microscopic and macroscopic one, give the same solution in term of the moments apart from the numerical errors. In particular, the DSMC method gives solutions which oscillates around the exact solution of the problem. This implies that, the above renormalization of the mass only force the density to be closer to the value furnished by the exact solution of the problem. In fact, the moments equations furnish solutions which contain less fluctuations with respect to the DSMC method. The reason for this lower level of noise of the macroscopic equations is that only the perturbation from the equilibrium g is computed statistically and not the full solution as in the original DSMC method.
We discuss now the matching of the momentum and of the energy, which are done after the matching of the density. In the Monte Carlo setting these moments can be obtained by the following piecewise constant reconstructions where u n j is a vector representing the mean velocity in the three spatial directions. To match mean velocity and energy with those of equations (44-45) we then apply the transformation described in [6] which permits to get a new set of velocities V * i defined by which gives After the matching at time n, the next point in the time marching procedure described before, regards the collisions. This step has been already described in section 3.1. It uses equation (35) to compute from f * , which is now the distribution function value after the transport and the matching, the new distribution f n .
We now discuss the sixth point of the method: the computation of g n from f n . The perturbation g n is given by where once again we used the first order in time integrating factor scheme (35). The above expression tells us that the moments of g n j can be obtained as a contribution of two terms. The first term is obtained by reconstructing from the particles the moments as in (49, 51). The second term is obtained by integrating over the velocity space the analytic expression of the Maxwellian distribution. This implies directly that the moments related to the second part does not contain any statistical error. Finally, observe that in the limit ε → 0 the contribution of the perturbation g goes to zero. As a consequence, the method proposed resolves a macroscopic system which is as closer as the Knudsen number diminishes to the compressible Euler equation. In the limit of ε = 0, the scheme exactly solve the compressible Euler equations without any source of statistical error. Finally, the discretized moments equations (45) can be rewritten as As ∆t/ε grows, which means that the system approaches the equilibrium, the contribution of the kinetic term vanishes even though it is evaluated through particles. Observe that, this does not happen if we just compute the kinetic term ∇ x · vmg from the particles without considering the structure of the distribution function f . This dramatically decreases fluctuations when the Knudsen number is small. As a conclusion for this section, we discuss some different reconstructions of the moments starting from the particles. Instead of using the piecewise constant reconstructions (49-51), one can think to smoother reconstructions which can probably further diminish the statistical noise. This procedure will be alternative to the moving average algorithm used for computing g from g (43). In other words, we first smooth the contribution of the particles in the computation of the macroscopic quantities and we then after just apply the finite volume method to set of resulting moments equations without using average algorithms.
A general way to compute the integral over the velocity space is to sum over the particles with the use of weight functions. In this strategy, the value of the perturbation g will be given by where B ≥ 0 is a suitable weight function s.t.

Numerical results
In this section we report some numerical results for the method proposed on different test cases. First, we consider two shock problems and then, we perform an accuracy test using a smooth initial data with periodic boundary conditions. For the two shock tests, we compare the moment guided (MG) solution with a Monte Carlo (MC) solution which employs the first order exponential Runge-Kutta method together with a splitting technique between the transport and the collision parts. We report also in each figure the results for the limit compressible Euler equation using the same second order MUSCL scheme which has been used for the equilibrium part of the Moment Guided method. Finally, the reference solutions for the two shock test problems are obtained by using the same Monte Carlo method described above with the same number of cells, in which, however, the number of particles is such that the statistical noise is very small. All the tests considered are relative to unsteady problems. We did this choice with the scope of highlighting the fact that the method proposed in this paper is specifically designed for situations in which the classical variance reduction techniques which employs time averaging, which is the typical case of DSMC methods, cannot be used or turns out to be useless, since time-averaging leads to the same computational effort of using more particles in one single simulation.

Unsteady shock test
The first problem concerns the study of an unsteady shock. The figures 1 to 3 consider the same initial data for the density ̺ = 1, the mean velocity u = −1 and the temperature T = 1 for   The three figures 1-3 show a large reduction of fluctuations for all cases analyzed. Observe that the reduction of the fluctuations depends on the Knudsen number value. In particular in the case ε = 10 −4 the solution is very close to the limit solution. Observe also that the solution furnished by the Moment Guided method when ε = 10 −4 is closer to the shock with respect to the reference Monte Carlo solution. This can be interpreted as an over relaxation problem of the Moment Guided method. However, an in deep simulation analysis we did not report, shows that, the difference between the MC and the MG methods when ε = 10 −4 is due to the more diffusive behavior of the Monte Carlo method. This is caused by the large time steps ∆t allowed by the exponential method. These large time steps introduce a numerical diffusion in the schemes which is more important for the Monte Carlo method and less important for the Moment Guided one.

Sod shock tube
The second problem analyzed is the Sod tube test. Again this choice relies on the fact that the method is specifically aimed to the study of unsteady problems. In this case, we want to study the capability of the method in describing different waves with respect to the simple shock wave. The figures 4 to 7 consider the same initial data for the density ̺ = 1 upstream and ̺ = 0.125 downstream of the initial shock, the mean velocity u = 0 in all the domain and the temperature: T = 5 upstream and T = 4 downstream of the shock. Different initial Knudsen number values are considered in the test, which range from ε = 10 −2 to ε = 10 −4 while the number of cell is 200. In figures 4, 6 and 7 200 particles per cell are also used on average, more precisely 200 2 particles are present at the beginning of the simulation and then distributed accordingly to the density in each cell. The final time is T f in = 0.08, while the time step as previously is given by the minimum between the the ratio of ∆x over the maximum velocity of the particles and the ratio of ∆x over the larger eigenvalue of the compressible Euler equations. Again the Knudsen number value does not constrain the time step, thanks to the exponential method employed in the resolution of the collision integral.
As for the previous test each figure depicts the density, the mean velocity and the temperature from top to bottom, with the Monte Carlo solver (left) and the Moment Guided method (right). The reference solution is obtained through the Monte Carlo method in which the number of particles is taken equal to 2 10 7 with a solution averaged 5 times, this gives a final number of 10 8 particles employed for constructing the reference solution. The figures show good results for all ranges of Knudsen numbers in terms of reduction of fluctuations. In particular, we clearly see a strong reduction of fluctuations for ε = 10 −3 (figure 6), while for ε = 10 −4 (figure 7) the fluctuations are almost completely disappeared. In the case of ε = 10 −2 (figure 4), we see that even if the statistical error is smaller for the Moment Guided method than for the Monte Carlo method, both solutions are still far from the converged solution indicated by the dash-dotted line. Thus, we reported the same simulation in figure 5 in which 1000 particle per cell on average are employed for both methods. This figure permits to show the faster convergence of the Moment Guided method towards the reference solution also in the case of larger Knudsen. Finally, thanks to the high order discretization of the moments equations, observe that when the gas is close to the fluid limit ( figure 7), the shock is very well represented.

Accuracy test
In this final part, we report the results of a an error analysis with respect to the number of particles.
In this case the reference solution is obtained by averaging M independent realizations for the for the Monte Carlo method (subscript M C) and for the Moment Guided method (subscript M G) Two different reference solutions have been used since the two schemes present different discretization errors and thus we expect them to converge, when the number of particles goes to infinity, to different discretized (in time and space) solutions. The difference between the two limit solutions (Monte Carlo and Moment Guided Monte Carlo) is mainly due to the different numerical diffusion introduced by the two methods. This has been made clear for instance in the unsteady shock test case when ε = 10 −4 . In this case, as already explained, the Monte Carlo solution is more diffusive with respect to the Moment Guided one and thus, in order to measure only the statistical error, we cannot consider the same reference solution for the two schemes. The two reference solutions have been obtained by fixing the time step and mesh size and letting the number of particles grows to eliminate the statistical error. The number of samples used to compute the reference solution is on average 10 5 per cell, while the number of realizations is M=5. The number of spatial cells is 100, the time step is fixed to 10 −3 for all the values of the Knudsen number. Both reference solutions still contain space and time discretization errors but the amount of such errors does not change when the number of particles varies. Therefore, the error we measure by fixing the space and time discretization is only due to the stochastic nature of the method. In practice, in our test we measured the quantity where U j represents the reference solutions, N = 100 and M = 10. The following initial data have been used ̺(x, 0) = 1 + a ̺ sin 2πx L u(x, 0) = 1.5 + a u sin 2πx L (60) E(x, 0) = 2.5 + a W sin 2πx L The boundary conditions are periodic while the time interval in which the simulation has been performed is t ∈ [0, 5 × 10 −2 ]. The results of this test in log-log scale are shown in Figure 8. On the left, we reported the error for the Monte Carlo while on the right for the Moment Guided method. For the Monte Carlo method, as expected the stochastic error does not change with respect to the Knudsen number. At variance, for the Moment Guided method, the error decreases as the Knudsen number diminishes. In particular observe that for the test studied the stochastic error is between 8 and 12 times smaller of the error of Monte Carlo method for ε = 10 −2 and ε = 10 −3 while is 10 −4 smaller when ε = 10 −4 .

Conclusions
We have extended a new class of hybrid methods, first developed in [11], which aim at reducing the variance in Monte Carlo schemes to the general case of the Boltzmann equation. The key idea consists in closing the set of macroscopic moments equations through a coupling with the kinetic equation solved by means of Monte Carlo methods. In order to correctly close the macroscopic system, we need to constrain particle positions and velocities in such a way that the moments given by the solution of the Boltzmann equation exactly match the moments given by the solution of the macroscopic set of moment equations. The new Moment Guided method is constructed in such a way that it avoids the problem of stiff regimes or equivalently of very small Knudsen numbers. In particular, the method benefits of the asymptotic preserving property which permits to circumvent the small scale resolution capturing the limit solution of the problem. A third important achievement consists in progressively reduce the contribution of the kinetic solution when the Knudsn number tends to zero and to automatically obtain in the fluid limit an high order method for the compressible Euler equations without any stochastic contribution. The numerical results performed show large reductions of fluctuations in all regimes, from dense to rarefied, compared to Monte Carlo methods for the Boltzman equation. The reduction becomes stronger as we approach equilibrium. In addition when ε = 0 the stochastic error becomes automatically zero. The numerical convergence tests show that we effectively get better performances for the Moment Guided method with respect to pure Monte Carlo schemes for all the cases studied. In particular, for unsteady problems, in which time averaging techniques lose their efficiency the method developed in this paper seems very promising. The computational cost of the method is comparable to the cost of a traditional Monte Carlo solver with the addition of the cost of a macroscopic solver for the compressible Euler equations. This latter is usually computationally much less expensive than any type of solver applied to kinetic equations.
Currently, we are working on extensions of the present method using higher order closures of the moment hierarchy in order to solve a larger set of hydrodynamics equations and thus to additionally reduce the stochastic errors. We are also working to the extension of the method to plasma physics problems such as the Vlasov and the collisional Vlasov models.